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ABSTRACT 

This  article  presents  an  introduction  to  multiscale  and  stabilized  methods,  which  represent  unified 
approaches  to  modeling  and  numerical  solution  of  fluid  dynamic  phenomena.  Finite  element 
applications  are  emphasized  but  the  ideas  are  general  and  apply  to  other  numerical  methods  as  well. 
(They  have  been  used  in  the  development  of  finite  difference,  finite  volume,  and  spectral  methods, 
in  addition  to  finite  element  methods.)  The  analytical  ideas  are  first  illustrated  for  time-harmonic 
wave-propagation  problems  in  unbounded  fluid  domains  governed  by  the  Helmholtz  equation.  This 
leads  to  the  well-known  Dirichlet-to-Neumann  formulation.  A  general  treatment  of  the  variational 
multiscale  method  in  the  context  of  an  abstract  Dirichlet  problem  is  then  presented  which  is  applicable 
to  advective-diffusive  processes  and  other  processes  of  physical  interest,  ft  is  shown  how  the  exact 
theory  represents  a  paradigm  for  subgrid-scale  models  and  a  posteriori  error  estimation.  Hierarchical 
p-methods  and  bubble  function  methods  are  examined  in  order  to  understand  and,  ultimately, 
approximate  the  “fine-scale  Green’s  function”  which  appears  in  the  theory.  Relationships  among  so- 
called  residual-free  bubbles,  element  Green’s  functions,  and  stabilized  methods  are  exhibited.  These 
ideas  are  then  generalized  to  a  class  of  non-symmetric,  linear  evolution  operators  formulated  in  space- 
time.  The  variational  multiscale  method  also  provides  guidelines  and  inspiration  for  the  development 
of  stabilized  methods  (e.g.,  SUPG,  GLS,  etc.)  which  have  attracted  considerable  interest  and  have 
been  extensively  utilized  in  engineering  and  the  physical  sciences.  An  overview  of  stabilized  methods 
for  advective-diffusive  equations  is  presented.  A  variational  multiscale  treatment  of  incompressible 
viscous  flows,  including  turbulence  is  also  described.  This  represents  an  alternative  formulation  of 
Large  Eddy  Simulation  which  provides  a  simplified  theoretical  framework  of  LES  with  potential  for 
improved  modeling. 
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1.  Introduction 

Stabilized  methods  were  originally  developed  about  25  years  ago  and  reported  on  in  a  series 
of  conference  papers  and  book  chapters.  The  first  archival  journal  article  appeared  in  1982 
(Brooks  and  Hughes,  1982).  This  work  summarized  developments  up  to  1982  and  brought  to 
prominence  the  SUPG  formulation  (i.e.,  Streamline  Upwind  Petrov-Galerkin) .  It  was  argued 
that  stability  and  accuracy  were  combined  in  this  approach,  and  thus  it  represented  an 
improvement  over  classical  upwind,  artificial  viscosity,  central-difference,  and  Galerkin  finite 
element  methods.  Mathematical  corroboration  came  shortly  thereafter  in  the  work  of  Johnson, 
Navert  and  Pitkaranta  (1984).  Subsequently,  many  works  appeared  dealing  with  fundamental 
mathematical  theory  and  diverse  applications.  A  very  large  literature  on  stabilized  methods 
has  accumulated  in  the  process. 

In  1995  it  was  shown  by  Hughes  (1995)  that  stabilized  methods  could  be  derived  from 
a  variational  multiscale  formulation.  Subsequently,  the  multiscale  foundations  of  stabilized 
methods  have  become  a  focal  point  of  research  activities  and  have  led  to  considerable 
conceptual  and  practical  progress.  The  view  taken  in  this  work  is  that  the  basis  of  residual- 
based,  or  consistent,  stabilized  methods  is  a  variational  multiscale  analysis  of  the  partial 
differential  equations  under  consideration.  This  approach  combines  ideas  of  physical  modeling 
with  numerical  approximation  in  a  unified  way.  To  provide  motivation  for  the  developments 
which  follow  a  relevant  physical  example  will  be  described  first. 

Considerations  of  environmental  acoustics  are  very  important  in  the  design  of  high-speed 
trains.  In  Japan,  environmental  laws  limit  the  sound  pressure  levels  50  meters  from  the  tracks. 
The  Shinkansen,  or  “bullet  trains”,  obtain  electric  power  from  pantographs  in  contact  with 
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Figure  1.  Shroud  surrounding  a  pantograph. 


overhead  lines.  In  order  to  reduce  aerodynamic  loads  on  pantographs,  “shrouds”  have  been 
designed  to  deflect  airflow.  Photographs  of  pantographs  and  shrouds  are  shown  in  Figures  1  and 
2.  The  shroud  reduces  structural  loads  on  the  pantograph  and,  concomitantly,  reduces  acoustic 
radiation  from  the  pantograph,  but  it  also  generates  considerable  noise  in  the  audible  range. 
Research  studies  have  been  performed  to  determine  the  acoustic  signatures  of  a  shroud  design 
(see,  e.g.,  Holmes  et  al.,  1997).  In  Holmes  et  al.  (1997)  a  Large  Eddy  Simulation  is  performed 
to  determine  the  turbulent  fluid  flow  in  the  vicinity  of  a  shroud.  See  Figure  3  for  a  schematic 
illustration.  There  are  several  important  length  scales  in  such  a  calculation.  Among  them  are 
L,  the  characteristic  length  scale  of  the  domain  in  which  there  are  strong  flow  gradients  and 
turbulence;  h,  the  characteristic  length  scale  of  the  mesh  used  in  the  numerical  analysis;  and  l , 
the  characteristic  length  of  the  smallest  turbulent  eddy.  These  scales  are  widely  separated,  that 
is,  L  3>  h  3>  /,  emphasizing  the  importance  of  multiscale  phenomena.  The  results  of  a  fluid 
dynamics  calculation  are  presented  in  Figure  4.  Note  the  smooth,  braided  vortical  structure  in 
the  boundary  layer,  and  the  turbulence  inside  and  above  the  shroud.  Eddies  impinge  on  the 
downwind  face  of  the  shroud  and  roof  of  the  train  and  give  rise  to  significant  pressures,  and 
ultimately  to  considerable  noise  propagation.  A  Fourier  transform,  with  respect  to  time,  of  the 
pressure  field  on  the  surface  of  the  shroud  and  roof  is  shown  in  Figure  5.  Note  the  spots  of 
intense  pressure.  These  locations  fluctuate  as  functions  of  frequency.  In  order  to  determine  the 
radiated  pressure  in  the  far  field,  the  Fourier  transformed  fluid  flow  is  used  to  generate  the  so- 


Encyclopedia  of  Computational  Mechanics.  Edited  by  Erwin  Stein,  Rene  de  Borst  and  Thomas  J.R.  Hughes. 
©  2004  John  Wiley  &  Sons,  Ltd. 


5 


Figure  2.  Pantographs  in  withdrawn  and  deployed  configurations. 
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Figure  3.  Turbulent  flow  in  a  domain  surrounding  a  pantograph  shroud  on  the  roof  of  a  high-speed 

train. 


called  Lighthill  turbulence  tensor  (Lighthill,  1952,1954),  from  which  sources  can  be  determined. 
These  are  used  to  drive  the  acoustic  field,  which  is  determined  by  solving  the  Helmholtz 
equation  (i.e.,  the  time-harmonic  wave  equation).  A  boundary- value  problem  needs  to  be  solved 
for  each  frequency  of  interest  in  order  to  construct  the  sound  pressure  level  spectrum.  These 
problems  are  classified  as  “exterior  problems,”  involving  domains  of  infinite  extent.  In  order 
to  use  a  domain-based  numerical  procedure,  such  as  finite  elements,  an  artificial  boundary 
is  introduced  which  surrounds  the  region  containing  the  acoustic  sources.  See  Figure  6  for  a 
schematic  illustration.  The  solution  of  the  problem  posed  within  the  artificial  boundary  needs 
to  approximate  the  solution  of  original  infinite-domain  problem.  This  necessitates  inclusion 
of  a  special  boundary  condition  on  the  artificial  boundary,  in  order  to  transmit  outgoing 
waves  without  reflection.  Various  schemes  have  been  proposed.  The  characteristic  length-scale 
induced  by  the  artificial  boundary,  R,  is  of  the  order  of  L  in  the  most  effective  approaches.  The 
distance  to  a  point  of  interest  in  the  far-held,  D ,  is  usually  much  larger  than  R.  The  solution  at 
D  can  be  determined  by  the  solution  on  the  artificial  boundary,  which  is  determined  from  the 
near-held  numerical  solution.  The  length  scales,  oo,  induce  additional  multiscale 

considerations.  A  sound  pressure  level  spectrum  at  a  microphone  location  (i.e.,  D)  is  compared 
with  numerical  results  in  Figure  7.  For  detailed  description  of  procedures  used  for  aeroacoustic 
and  hydroacoustic  applications,  see  Oberai,  Roknaldin  and  Hughes  (2000,2002). 

The  multiscale  aspects  of  the  fluid  how  and  acoustic  propagation  will  be  discussed  in 
the  sequel.  The  former  problem  is  nonlinear  and  more  complex.  It  will  be  described  in  the 
last  section.  In  the  next  section,  a  multiscale  formulation  of  acoustic  radiation  is  presented. 
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Figure  4.  Turbulent  flow  about  a  pantograph  shroud. 


The  connections  between  multiscale  formulations  and  stabilized  methods  are  developed  in  the 
intervening  sections. 

The  terminology  “multiscale”  is  used  widely  for  many  different  things.  Other  concepts  of 
multiscale  analysis  are,  for  example,  contained  in  E  and  Engquist  (2003)  and  Wagner  and  Liu 
(2003). 


2.  Dirichlet-to-Neumann  Formulation 

The  exterior  problem  for  the  Helmholtz  equation  (i.e.,  the  complex- valued,  time-harmonic, 
wave  equation)  is  considered.  The  viewpoint  adopted  is  that  there  are  two  sets  of  scales  present, 
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Figure  5. 


one  associated  with  the  near  field  and  one  associated  with  the  far  field.  The  near-field  scales 
are  viewed  as  those  of  the  exact  solution  exterior  to  the  body,  but  within  an  enclosing  simple 
surface,  such  as,  for  example,  a  sphere.  The  enclosing  surface  is  not  part  of  the  specification  of 
the  boundary-value  problem,  but  rather  it  is  specified  by  the  analyst.  The  near-field  scales  are 
viewed  as  numerically  “resolvable”  in  this  case.  They  may  also  be  thought  of  as  local  or  small 
scales.  The  scales  associated  with  the  solution  exterior  to  the  sphere  (far  field)  are  the  global  or 
large  scales  and  are  viewed  as  numerically  “unresolvable”  in  the  sense  that  the  infinite  domain 
of  the  far  field  cannot  be  dealt  with  by  conventional  bounded-domain  discretization  methods. 
The  solution  of  the  original  problem  is  decomposed  into  non-overlapping  near-field  and  far-held 
components,  and  the  far-held  component  is  exactly  solved  for  in  terms  of  the  exterior  Green’s 
function  satisfying  homogeneous  Dirichlet  boundary  conditions  on  the  sphere.  (Shapes  other 
than  a  sphere  are  admissible,  and  useful  in  particular  cases,  but  for  each  shape  one  must 
be  able  to  solve  the  exterior  Green’s  function  problem  in  order  to  determine  the  far-held 
solution.)  The  far-held  component  of  the  solution  is  then  eliminated  from  the  problem  for  the 
near  held.  This  results  in  a  well-known  variational  formulation  on  a  bounded  domain  which 
exactly  characterizes  the  near-held  component  of  the  original  problem.  It  is  referred  to  as  the 
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Figure  6.  A  numerical  problem  is  solved  within  the  artificial  boundary  at  radius  R.  The  point  of 
interest  is  located  at  D.  The  analytical  problem  involves  a  domain  of  infinite  extent. 


Figure  7.  Sound  pressure  level  spectrum. 
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n 


Figure  8.  An  exterior  domain. 


Dirichlet-to-Neumann  (DtN)  formulation  because  of  the  form  of  the  boundary  condition  on 
the  sphere  in  the  problem  on  the  bounded  domain  (Givoli,  1992;  Givoli  and  Keller,  1988,1989; 
Harari  and  Hughes,  1992,1994).  The  so-called  DtN  boundary  condition  is  nonlocal  in  the  sense 
that  it  involves  an  integral  operator  coupling  all  points  on  the  sphere.  Nonlocality  is  a  typical 
ingredient  in  formulations  of  multiscale  phenomena. 


2.1.  Dirichlet-to-Neumann  formulation  for  the  Helmholtz  operator 

Consider  the  exterior  problem  for  the  Helmholtz  operator.  Let  Q  C  be  an  exterior  domain, 
where  d  is  the  number  of  space  dimensions  (see  Figure  8).  The  boundary  of  O  is  denoted  by  T 
and  admits  the  decomposition 


r  =  rg\Jrh  (l) 

0  =  r9p|r„  (2) 

where  Tg  and  T ft  are  subsets  of  T.  The  unit  outward  vector  to  T  is  denoted  by  n.  The  boundary- 
value  problem  consists  of  finding  a  function  u  :  — >  C,  such  that  for  given  functions  /  :  O  — >  C, 
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n 


Figure  9.  Decomposition  of  12  into  a  bounded  domain  12  and  an  exterior  domain  12'. 


g  :  rg  — y  C  and  h  :  Th  — »  C,  the  following  equations  are  satisfied: 


Cu  =  f 
■u  =  g 

u,n  =  ikh 

d—1 

lim  r  2  (u  r  —  iku)  =  0 

r— ►  oo 

where 

-£  =  A  +  A:2 


in  0 

(3) 

on  Tg 

(4) 

on  r  h 

(5) 

(Sommerfeld  radiation  condition) 

(6) 

(Helmholtz  operator) 

(7) 

and  k  G  C  is  the  wave  number,  i  =  y/—T,  and  A  is  the  Laplacian  operator.  The  radial 
coordinate  is  denoted  by  r  and  a  comma  denotes  partial  differentiation.  The  Sommerfeld 
radiation  condition  enforces  the  condition  that  waves  at  infinity  are  outgoing. 

Next,  consider  a  decomposition  of  the  domain  O  into  a  bounded  domain  O  and  an  exterior 
region  O'.  The  boundary  which  separates  O  and  O'  is  denoted  T^.  It  is  assumed  to  have  a 
simple  shape  (e.g.,  spherical).  See  Figure  9.  The  decomposition  of  O,  and  a  corresponding 
decomposition  of  the  solution  of  the  boundary-value  problem,  are  expressed  analytically  as 
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follows: 

0  = 

n  (J  iv 

(8) 

0  = 

np|sr 

(9) 

u  = 

u  +  v! 

(sum  decomposition) 

(10) 

U  In'  = 

u'  In 

o  o 

(disjoint  decomposition) 

(11) 

u  = 

{:■ 

on  f \ 
on  f Y 

(12) 

Think  of  u  as  the  near-field  solution  and  v!  as  the  far-held  solution. 


2.2.  Exterior  Dirichlet  problem  for  u' 

Attention  is  now  focused  on  the  problem  in  the  domain  O'  exterior  to  T#.  The  unit  outward 
normal  vector  on  T^  (with  respect  to  f V)  is  denoted  n'  (see  Fig.  10).  Assume  that  /  vanishes 
in  the  far  held,  that  is, 

/  =  0  on  it!  (13) 

The  exterior  Dirichlet  problem  consists  of  finding  a  function  v!  :  O  — >  C  such  that 

Eu!  =  0  in  0!  (14) 

u!  =  u  onT#  (15) 

lim  r-2""  {y!  —  iku')  =  0  (16) 

r — >-oo 

Note  that  the  boundary  condition  (15)  follows  from  the  continuity  of  u  across  r r. 

2.3.  Green’s  function  for  the  exterior  Dirichlet  problem 

The  solution  of  the  exterior  Dirichlet  problem  can  be  expressed  in  terms  of  a  Green’s  function 
g  satisfying 


Cg  =  5  in  Q' 

(17) 

g  =  0  on  rR 

(18) 

From  Green’s  identity, 

d-1 

lim  r  2  (g  r  —  ikg )  =  0 

1 — >oo 

t* 

(19) 

u'(.y)  =  -  g  ,n'x{x,y)u'(x)dTx 

JrR 

(20) 

The  so-called  DtN  map  is 

obtained  from  (20)  by  differentiation  with  respect  to  n', 

u\n'(y)  = 

-  g  ,n'n'(x,y)u'(x)dTx  aa  u',(y)  =  Mu' 

JrR 

(21) 
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Figure  10.  Domain  for  the  far-held  problem. 


The  DtN  map  is  used  to  develop  a  formulation  for  u  on  the  bounded  domain  fi.  In  this  way 
the  far-held  phenomena  are  incorporated  in  the  problem  for  the  near  held.  The  Dirichlet-to- 
Neumann  formulation  for  u  will  be  developed  by  way  of  a  variational  argument. 

Let  II  denote  the  potential  energy  for  the  original  boundary- value  problem,  namely 


where 


II(u)  =  n(u  +  u') 

=  a(u ,  u)  +  ^ a(u v!)  -  (u,  /)  -  (u,  ikh) r 


a(w,  u ) 

K  /) 

(w,  ikh) r 


/  (Vto  •  Vu  -  k2wu)dn 
Jn 


wfdfl 


w  ikh  dT 


lrh 


Consider  a  one-parameter  family  of  variations  of  u,  that  is, 

(u  +  u')  +  e(w  +  w') 


subject  to  the  following  continuity  constraints 


u  =  v!  on  T/j 

w  =  w'  on  Tfl 


(22) 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


Encyclopedia  of  Computational  Mechanics.  Edited  by  Erwin  Stein,  Rene  de  Borst  and  Thomas  J.R.  Hughes. 
©  2004  John  Wiley  &  Sons,  Ltd. 


14 


ENCYCLOPEDIA  OF  COMPUTATIONAL  MECHANICS 


where  e  G  I  is  a  parameter.  Taking  the  Frechet  derivative,  the  first  variation  of  II  is  calculated 
as  follows: 

0  =  DU{u  +  u')  ■  (w  +  w') 

=  a(w,  u)  +  a(w',  u)  —  (w,  /)  —  ( w ,  ikh)  r 

=  a(w,u)  +  (wr ,Cu')  +  (w',u'n,)rR  —  (w,  /)  —  ( w,ikh)r 

=  a(w,  u)  +  0  +  (w,  Mu)yr  —  (w,  f)  —  (w,  ikh)r  (29) 

where 

(w,Mu)Vr=  1  /  w(y)gtnxnv(x,y)  u(x)dYxdTy  (30) 

Jtr  Jtr 

In  obtaining  (29),  (21)  and  the  continuity  conditions,  (27)  and  (28),  have  been  used.  Note 
that  in  (30)  differentiation  with  respect  to  n  =  —n'  has  been  employed.  Equation  (29)  can  be 
written  concisely  as 

(31) 

where 

B(w,u;  g)  =  a(w,u)  +  (w,  Mu)rR  (32) 

L(w)  =  (w,f)  +  {w,ikh)r  (33) 

Remarks 

1.  (31)  is  an  exact  characterization  of  u. 

2.  The  effect  of  u!  on  the  problem  for  u  is  nonlocal.  The  additional  term,  (30),  is  referred  to 
as  the  DtN  boundary  condition.  It  represents  a  perfect  interface  that  transmits  outgoing 
waves  without  reflection. 

3.  (31)  is  the  basis  of  numerical  approximations,  viz. 

B(wh,uh;g)  =  L(Wh)  (34) 

where  wh  and  uh  are  finite-dimensional  approximations  of  w  and  u,  respectively. 

4.  In  practice,  M  (or  equivalently  g)  is  also  approximated  by  way  of  truncated  series, 
differential  operators,  etc.  Thus,  in  practice,  we  work  with 

B(wh,uh;g)  =  L(wh)  (35) 

where 

g  «  g  (36) 

2-4-  Bounded  domain  problem  for  u 

The  Euler- Lagrange  equations  of  the  variational  formulation  give  rise  to  the  boundary- 
value  problem  for  u  on  the  bounded  domain  f l  (see  Fig.  11),  that  is, 


Cm 

f 

in  O 

(37) 

si  = 

9 

on  rg 

(38) 

U,n  = 

ikh 

on 

(39) 

U,n  = 

—Mu 

on  Tr 

(40) 
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Figure  11.  Bounded  domain  for  the  near-field  problem. 


The  preceding  developments  may  be  summarized  in  the  following  statements: 

1.  u  =  u  +  v!  ( disjoint  sum  decomposition). 

2.  v!  is  determined  analytically. 

3.  u'  is  eliminated,  resulting  in  a  formulation  for  u  which  is  the  basis  of  numerical 
approximations . 

4.  The  effect  of  u'  is  nonlocal  in  the  problem  for  u. 

5.  Interpreted  as  a  multiscale  problem,  u'  represents  the  large  scales  of  the  far  field, 
whereas  u  represents  the  small  scales  of  the  near  field. 


3.  Variational  Multiscale  Method 

The  variational  multiscale  method  is  a  procedure  for  deriving  models  and  numerical  methods 
capable  of  dealing  with  multiscale  phenomena  ubiquitous  in  science  and  engineering.  It  is 
motivated  by  the  simple  fact  that  straightforward  application  of  Galerkin’s  method  employing 
standard  bases,  such  as  Fourier  series  and  finite  elements,  is  not  a  robust  approach  in  the 
presence  of  multiscale  phenomena.  The  variational  multiscale  method  seeks  to  rectify  this 
situation.  The  anatomy  of  the  method  is  simple:  sum  decompositions  of  the  solution, 
u  =  u+u' ,  are  considered  where  u  is  solved  for  numerically .  An  attempt  is  made  to  determine 
u!  analytically ,  eliminating  it  from  the  problem  for  u.  u  and  v!  may  overlap  or  be  disjoint , 
and  u'  may  be  globally  or  locally  defined.  The  effect  of  u'  on  the  problem  for  u  will  always 
be  nonlocal.  In  the  previous  section,  the  variational  multiscale  method  was  used  to  derive 
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the  Dirichlet-to-Neumann  formulation  of  the  Helmholtz  equation  in  an  unbounded  domain.  In 
this  section,  attention  is  confined  to  cases  on  bounded  domains  in  which  u  represents  “coarse 
scales”  and  u'  “fine  scales”. 

An  attempt  is  made  to  present  the  big  picture  in  the  context  of  an  abstract  Dirichlet 
problem  involving  a  second-order  differential  operator  which  is  assumed  nonsymmetric  and/or 
indefinite.  This  allows  consideration  of  equations  of  practical  interest,  such  as  the  advection- 
diffusion  equation,  a  model  for  fluid  mechanics  phenomena,  and  the  Helmholtz  equation,  of 
importance  in  acoustics  and  electromagnetics.  After  introducing  the  variational  formulation  of 
the  Dirichlet  problem,  its  multiscale  version  is  described. 

First  the  “smooth  case”  is  considered,  in  which  it  is  assumed  that  all  functions  are  sufficiently 
smooth  so  that  distributional  effects  (e.g.,  Dirac  layers)  may  be  ignored.  This  enables  a  simple 
derivation  of  the  exact  equation  governing  the  coarse  scales.  It  is  helpful  to  think  of  this  case  as 
pertaining  to  the  situation  in  which  both  the  coarse  and  fine  scales  are  represented  by  Fourier 
series. 

Next,  a  case  of  greater  practical  interest  is  considered  in  which  standard  finite  elements  are 
employed  to  represent  the  coarse  scales.  Due  to  lack  of  continuity  of  derivatives  at  element 
interfaces,  it  is  necessary  to  explicitly  account  for  the  distributional  effects  omitted  in  the 
smooth  case.  This  is  referred  to  as  the  “rough  case”.  Again,  an  exact  equation  is  derived 
governing  the  behavior  of  coarse  scales.  It  is  this  equation  that  is  proposed  as  a  paradigm 
for  developing  subgrid-scale  models.  Two  distinguishing  features  characterize  this  result.  The 
first  is  that  the  method  may  be  viewed  as  the  classical  Galerkin  method  plus  an  additional 
term  driven  by  the  distributional  residual  of  the  coarse  scales.  This  involves  residuals  of  the 
partial  differential  equation  under  consideration  on  element  interiors  (this  is  the  smooth  part 
of  the  residual),  and  jump  terms  involving  the  boundary  operator  on  element  interfaces  (this 
is  the  rough  part  deriving  from  Dirac  layers  in  the  distributional  form  of  the  operator).  The 
appearance  of  element  residuals  and  jump  terms  are  suggestive  of  the  relationship  between 
the  multiscale  formulation  and  various  stabilized  methods  proposed  previously.  The  second 
distinguishing  feature  is  the  appearance  of  the  fine-scale  Green’s  function.  In  general,  this 
is  not  the  classical  Green’s  function,  but  one  that  emanates  from  the  fine-scale  subspace. 
It  is  important  to  note  that  the  fine-scale  subspace,  V' ,  is  infinite-dimensional,  but  a  proper 
subspace  of  the  space,  V,  in  which  it  is  attempted  to  solve  the  problem.  The  direct  sum 
relationship  V  =  V  ©  V'  where  V  is  the  coarse-scale,  finite  element  subspace  is  satisfied. 
A  problem  that  arises  in  developing  practical  approximations  is  that  the  fine-scale  Green’s 
function  is  nonlocal. 

Before  addressing  this  issue,  the  relationship  between  the  fine-scale  solution  and  a  posteriori 
error  estimation  is  discussed.  It  is  noted  first  that  by  virtue  of  the  formulation  being  exact, 
the  fine-scale  solution  is  precisely  the  error  in  the  coarse-scale  solution.  Consequently,  the 
representation  obtained  of  the  fine-scale  solution  in  terms  of  the  distributional  coarse-scale 
residual  and  the  fine-scale  Green’s  function  is  a  paradigm  for  a  posteriori  error  estimation.  It 
is  then  noted  that  it  is  typical  in  a  posteriori  error  estimation  procedures  to  involve  the  element 
residuals  and/or  interface  jump  terms  as  driving  mechanisms.  The  mode  of  distributing  these 
sources  of  error  may  thus  be  inferred  to  be  approximations  of  the  fine-scale  Green’s  function.  As 
a  result,  it  is  clear  that  in  a  posteriori  error  estimation,  the  proper  distribution  of  residual  errors 
strongly  depends  on  the  operator  under  consideration.  In  other  words,  there  is  no  universally 
appropriate  scheme  independent  of  the  operator.  (A  similar  observation  may  be  made  for 
subgrid-scale  models  by  virtue  of  the  form  of  the  coarse-scale  equation.)  The  implications  of 
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Variational  form  Multiscale  Exact  equation 

of  PDE  exhibiting  analysis  for  ,7 


Figure  12.  The  variational  multiscale  method  is  a  framework  for  the  construction  of  subgrid- 
scale  models  and  effective  numerical  methods  for  partial  differential  equations  exhibiting  multiscale 
phenomena.  It  provides  a  physical  context  for  understanding  methods  based  on  residual-free  bubbles 

and  stabilized  methods. 


the  formula  for  the  fine-scale  solution  with  respect  to  a  posteriori  error  estimation  for  finite 
element  approximations  of  the  advection-diffusion  and  Helmholtz  equations  are  discussed. 

Next,  hierarchical  p-refinement  and  bubbles  are  examined  in  an  effort  to  better  understand 
the  nature  of  the  fine-scale  Green’s  function  and  to  deduce  appropriate  forms.  V  is  identified 
with  standard,  low-order  finite  elements,  and  V  with  the  hierarchical  basis.  An  explicit 
formula  for  the  fine-scale  Green’s  function  in  terms  of  the  hierarchical  basis  is  derived.  It  is 
concluded  that,  despite  the  nonlocal  character  of  the  fine-scale  Green’s  function,  it  can  always 
be  represented  in  terms  of  a  finite  basis  of  functions  possessing  local  support.  In  one-dimension, 
this  basis  consists  solely  of  bubbles,  in  two  dimensions,  bubbles  and  edge  functions;  etc.  This 
reduces  the  problem  of  approximating  the  Green’s  function  to  one  of  obtaining  a  good-quality, 
finite-dimensional  fine-scale  basis.  This  becomes  a  fundamental  problem  in  the  construction 
of  practical  methods.  Once  solved,  a  subgrid-scale  model  governing  the  coarse-scales,  and  an 
approximate  representation  of  the  fine-scale  solution  which  does  double  duty  as  an  a  posteriori 
error  estimator  for  the  coarse-scale  solution,  are  obtained. 

What  constitutes  a  good-quality,  but  practical,  fine-scale  basis  is  described  by  reviewing  the 
concept  of  residual- free  bubbles  (see  Baiocchi,  Brezzi  and  Franca,  1993).  The  use  of  fine-scale 
Green’s  functions  supported  by  individual  elements  is  then  reviewed.  Residual-free  bubbles 
and  element  Green’s  functions  are  intimately  related  as  shown  in  Brezzi  et  al.  (1997).  These 
concepts  may  be  used  to  derive  stabilized  methods  and  identify  optimal  parameters  which 
appear  in  their  definition.  The  ideas  are  illustrated  with  one-dimensional  examples. 
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Figure  13.  Domain  and  boundary  for  the  abstract  Dirichlet  problem. 


This  section  is  concluded  with  a  summary  of  results  and  identification  of  some  outstanding 
issues.  The  overall  flow  of  the  main  relationships  is  presented  in  Figure  12. 

An  alternative  approach  for  constructing  a  fine-scale  basis  can  be  found  in  the  literature 
describing  the  Discontinuous  Enrichment  Method  (DEM)  with  Lagrange  multipliers,  Farhat, 
Harari  and  Franca  (2001),  Farhat,  Harari  and  Hetmaniuk  (2003a,  2003b),  and  Harari,  Farhat 
and  Hetmaniuk  (2003).  In  this  hybrid  variational  multiscale  approach,  the  fine-scales  are 
based  on  the  free-space  solutions  of  the  homogeneous  differential  equation  to  be  solved.  For 
example,  for  the  Helmholtz  equation,  these  scales  are  represented  analytically  by  plane  waves. 
This  approach  leads  to  fine-scales  that,  unlike  bubbles,  do  not  vanish  but  are  discontinuous 
on  the  element  boundaries.  This  allows  circumventing  both  the  difficulty  in  attempting  to 
approximate  the  global  fine-scale  Green’s  function,  and  the  loss  of  some  global  effects  due  to 
the  restriction  of  residual-free  bubbles  to  a  vanishing  trace  on  the  element  boundaries.  However, 
the  DEM  approach  for  constructing  a  fine-scale  basis  introduces  additional  unknowns  at  the 
element  interfaces  in  the  form  of  Lagrange  multipliers  to  enforce  a  weak  continuity  of  the 
solution. 

3.1.  Abstract  Dirichlet  problem 

Let  fi  C  Rd,  where  d  >  1  is  the  number  of  space  dimensions,  be  an  open  bounded  domain 
with  smooth  boundary  T  (see  Fig.  13).  Consider  the  following  boundary- value  problem:  find 
u  :  O  — ►  ffi  such  that 


Cu  =  f  in  fi  (41) 

u  =  g  on  T  (42) 

where  /  :  O  — >  R  and  g  :  T  — >  ffi  are  given  functions.  Think  of  £  as  a  second-order  and,  in 
general,  nonsymmetric  differential  operator. 

3.1.1.  Variational  formulation  Let  S  C  H1( Cl)  denote  the  trial  solution  space  and 
V  C  1L1(H)  denote  the  weighting  function  space ,  where  iL1(H)  is  the  Sobolev  space  of 
square-integrable  functions  with  square-integrable  derivatives.  Assume  that  S  and  V  possess 
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Figure  14.  Coarse  and  fine  scale  components. 


the  following  properties: 

u  =  g  on  T  Vit  £  S  (43) 

w  =  0  on  T  \/w  €  V  (44) 

The  variational  counterpart  of  the  boundary-value  problem  (41)-(42)  is  given  as  follows:  find 
u  £  S  such  that  \/w  £  V 

a(w,  u )  =  (w,  f)  (45) 

where  (•,  •)  is  the  L,2{CL)  inner  product,  and  a(-,  •)  is  a  bilinear  form  satisfying 

a(w,  u)  =  (w,  Cu)  (46) 

for  all  sufficiently  smooth  w  £  V  and  u  £  S. 

3.2.  Variational  multiscale  method 
Let 


u  =  u  +  u  ( overlapping  sum  decomposition)  (47) 

where  u  represents  coarse  scales  and  u'  represents  fine  scales  (see  Fig.  14).  Likewise,  let 

w  =  w  +  w' .  (48) 

Let  S  =  S  ®  S'  and  V  =  V  ®  V'  where  S  (resp.,  S')  is  the  trial  solution  space  for  coarse 
(resp.,  fine)  scales  and  V  (resp.,  V')  is  the  weighting  function  space  for  coarse  (resp.,  fine) 
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scales.  Assume 


u 

=  9 

on  r 

Vues 

(49) 

v! 

=  0 

on  r 

Vu’  e  S' 

(50) 

w 

=  0 

on  r 

Vw  e  V 

(51) 

w' 

=  0 

on  r 

Vw'  e  v' 

(52) 

We  assume  S'  =  V' .  The  objective  is  to  derive  an  equation  governing  u. 


Remarks 

1.  It  is  helpful  to  think  of  S  and  V  as  finite-dimensional,  whereas  S'  and  V'  are  necessarily 
infinite-dimensional. 

2.  In  order  to  make  the  notion  of  the  direct  sums  precise,  one  needs  to  introduce  projectors 
II5  :  S  — >  S  and  Ily  :  V  — »  V,  such  that  u  =  II5M,  w  =  Ilyin,  11^  =  id  —  II5, 
Ily  =  id  —  Ily,  and,  in  particular,  i Ifj  =  Ily  :=  11'. 

3.2.1.  Smooth  case  The  developments  are  begun  by  considering  the  case  in  which  all  functions 


are  smooth.  The  idea  for  u  =  u  +  u'  is  illustrated  in  Figure  15.  The  situation  for  w  =  w  +  w' 
is  similar.  Assume  the  following  integration-by-parts  formulas  hold: 

a(w,u')  =  (jC*w,u')  Vw  e  V,  u!  e  S'  (53) 

a(w',u)  =  ( w',Cu )  Vw'eV',ueS  (54) 

a(w',u')  =  ( w',Cu' )  Vw'  e  V' ,  u'  €  S'  (55) 

Exact  variational  equation  for  u  (smooth  case)  Substitute  (47)  and  (48)  into  (45): 

a(w  +  w',u  +  u')  =  (w  +  w',  f)  Vw  e  V,  Viz/  €  V  (56) 

By  virtue  of  the  linear  independence  of  w  and  w' ,  (56)  splits  into  two  problems: 

Problem  (1)  a(w,u)  +  a(w,u')  =  (w,  f)  \/wCV  (57) 

a(w,  u)  +  ( C*w ,  u)  =  (■ w ,  /)  (58) 

Problem  (2)  a(w',u)  +  a(w',u')  =  (iz/ ,  /)  Viz/  €  V  (59) 

(w' ,  Cu)  +  (w1 ,  Cu')  =  (w' ,  f)  (60) 

In  arriving  at  (58)  and  (60),  the  integration-by-parts  formulas  (53)-(55)  have  been  employed. 
Rewrite  (60)  as 

(n')*^'  =  -(n ')*(£u-f)  in  Cl  (61) 

u'  =  0  on  T  (62) 


where  (IT)*  denotes  projection  onto  (V/)*,  the  dual  space  of  V'.  Endeavor  to  solve  this  problem 
for  v!  and  eliminate  u'  from  the  equation  for  u ,  namely  (58).  This  can  be  accomplished  with 
the  aid  of  a  Green’s  function. 
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Figure  15.  The  case  in  which  u  and  u'  are  smooth. 


Green’s  function  for  the  “dual  problem”  Consider  the  following  Green’s  function 
problem  for  the  adjoint  operator: 

(n,)*£*n/g(s,i/)  =  (n')*6(z-y)  VxeQ  (63) 

g (x,y)  =0  Vcc  G  T  (64) 

we  seek  agi  ker(  (n,)*£*n' ).  Let  g'  =  n^LL)*.  In  terms  of  the  solution  of  this  problem,  u' 
can  be  expressed  as  follows: 

u'{y)  =  -  [  g'{x,  y)  ( Cu  -  /)  (a;)  dClx  (65) 

J  n 

Equivalently,  (65)  can  be  written  in  terms  of  an  integral  operator  M'  as 

u'  =  M'  {Cu  -  /)  (66) 


Remarks 

1.  Cu  —  f  is  the  residual  of  the  coarse  scales. 

2.  The  fine  scales,  v! ,  are  driven  by  the  residual  of  the  coarse  scales. 

3.  It  is  very  important  to  observe  that  g'  is  not  the  usual  Green’s  function  associated 
with  the  corresponding  strong  form  of  (63).  Rather,  g'  is  defined  entirely  in  terms  of 
the  space  of  fine  scales,  namely  V' ■  Later  on,  an  explicit  formula  for  g'  will  be  derived 
in  terms  of  a  basis  for  V'. 
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Substituting  (66)  into  (58)  yields 


a(w,  u)  +  ( C*w ,  M'  ( Cu  —  /))  =  ( w ,  /)  \/w  €  V' 


where,  from  (65), 

(C*w,  M'  (Cu  —  /))  =  —  f  f  (C*w)(y)g'(x,y)(Cu- f)(x)dClx  d£Ly 
Jn  Jn 

Remarks 


(67) 


(68) 


1.  This  is  an  exact  equation  for  the  coarse  scales. 

2.  The  effect  of  the  fine  scales  on  the  coarse  scales  is  nonlocal. 

3.  By  virtue  of  the  smoothness  assumptions,  this  result  is  appropriate  for  spectral  methods, 
or  methods  based  on  Fourier  series,  but  it  is  not  sufficiently  general  as  a  basis  for  finite 
element  methods.  In  what  follows,  the  smoothness  assumption  is  relaxed  and  the  form 
of  the  coarse-scale  equation  appropriate  for  finite  elements  is  considered. 

3.2.2.  Rough,  case  (FEM)  Consider  a  discretization  of  O  into  finite  elements.  The  domain 
and  boundary  of  element  e,  where  e  G  {1,  2,  •  •  •  ,  nei},  in  which  ne i  is  the  number  of  elements, 
are  denoted  Oe  and  Te,  respectively  (see  Fig.  16).  The  union  of  element  interiors  is  denoted  f V 
and  the  union  of  element  boundaries  modulo  T  (also  referred  to  as  the  element  interfaces 
or  skeleton)  is  denoted  T',  viz. 


O' 

r' 

o 


Un‘ 

e=l 

(riel  \ 

y,Tr 

closure  (O') 


(69) 

(70) 

(71) 


Let  S,Vc  C°(0)  fl  H1( O)  be  classical  finite  element  spaces.  Note  that  S'  =  V'  C  H1( O), 
but  is  otherwise  arbitrary.  In  this  case  u  and  w  are  smooth  on  element  interiors  but  have  slope 
discontinuities  across  element  boundaries  (see,  e.g.,  Fig.  17). 

It  is  necessary  to  introduce  some  terminology  used  in  the  developments  which  follow.  Let 
(•,-)w  be  the  L2(co)  inner  product  where  oj  =  O,  Oe,  Te,  O',  T',  etc.  Recall,  (•,•)  =  (-,-)f2-  Let  [•] 
denote  the  jump  operator ,  viz.,  if  v  is  a  vector  field  experiencing  a  discontinuity  across  an 
element  boundary  (e.g.,  v  =  'S/w1w  €  V),  then 


n  ■  u] 

=  7l+ 

•  v+  +  n~ 

•  V 

=  n+ 

■  v+  —  n+ 

■  V 

—  n  •  i 

(■u+  —  v~), 

(72) 


where 


n  =  n+  =  —  n 


(73) 
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Figure  16.  Discretization  of  11  into  element  subdomains. 


Figure  17.  u  is  the  piecewise  linear  interpolate  of  u. 


is  a  unit  normal  vector  on  the  element  boundary  and  the  ±  designations  are  defined 
as  illustrated  in  Figure  18.  Note  that  (72)  is  invariant  with  respect  to  interchange  of  ± 
designations. 

In  the  present  case  there  is  smoothness  only  on  element  interiors.  Consequently,  integration- 
by-parts  gives  rise  to  nonvanishing  element  boundary  terms.  For  example,  if  w  £  V  and  u'  €  S' , 
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Figure  18.  Definition  of  unit  normals  on  an  element  boundary. 


O  X 


V 

V 

Figure  19.  Generalized  derivatives  of  piecewise  linear  and  quadratic  finite  elements. 
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the  following  integration-by-parts  formula  holds 

He  1 

a(w,u')  =  ((C*w,  u')Qe  +  u')re) 

e=l 

=  (£*w,  +  ([&*«>],«% 

=  (£*w,u%  (74) 

where  b*  is  the  boundary  operator  corresponding  to  C*  (e.g.,  if  C*  =  C  =  —A,  then 
6*  =  b  =  n  ■  V  =  d/dn).  Note,  from  (74),  there  are  three  different  ways  to  express  the 
integration-by-parts  formula.  The  first  line  of  (74)  amounts  to  performing  integration-by-parts 
on  an  element-by-element  basis.  In  the  second  line,  the  sum  over  element  interiors  has  been 
represented  by  integration  over  Q,'  and  the  element  boundary  terms  have  been  combined  in 
pairs,  the  result  being  a  jump  term  integrated  over  element  interfaces.  Finally,  in  the  third 
line,  C*w  is  viewed  as  a  Dirac  distribution  defined  on  the  entire  domain  O.  To  understand 
this  interpretation,  consider  the  following  example: 

Let 

C*w  =  wtXX  (75) 

and  assume  w  consists  of  piecewise  linear,  or  quadratic,  finite  elements  in  one  dimension.  The 
set-up  is  illustrated  in  Figure  19.  Note  that  w)XX  consists  of  Dirac  delta  functions  at  element 
boundaries  and  smooth  functions  on  element  interiors.  This  amounts  to  the  distributional 
interpretation  of  C*w  in  the  general  case.  It  is  smooth  on  element  interiors  but  contains  Dirac 
layers  on  the  element  interfaces,  which  give  rise  to  the  jump  terms  in  the  second  line  of  (74). 

Likewise,  there  are  additional  integration-by-parts  formulas:  \/w'  £  V',  u  £  S  and  u'  £  S' , 

"el 

a(w',u)  =  ((w',  Cu)ne  +  (w',  bu)re) 

e—l 

=  (w',Cu)n,  +  (if/,  [buj)r, 

=  (w',£u)n  (76) 


riel 

a(w',u')  =  Y,  (( w Cu')Qe  +  (w',  bu')re ) 

e—l 

=  (w',Cu')Q,  +  («/,  [Zm'J)r, 

=  (w/,  Cu')n  (77) 

where,  again,  £.u  and  Cu'  are  Dirac  distributions  on  fi. 

Exact  variational  equation  for  u  (rough  case)  The  distributional  interpretation  of  Cu, 
Cu'  and  C*w  allows  one  to  follow  the  developments  of  the  smooth  case  (see  Section  3.2.1). 
Thus,  the  formula  for  u'  can  be  expressed  in  three  alternative  forms  analogous  to  those  of  the 
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integration-by-parts  formulas,  viz., 

u'{y )  =  -  [  g '{x,y){Cu- f){x)dttx 

Jn 

=  -  f  g '{x,y){Cu~  f){x)dClx  -  f  g\x,y){bu\{x)dTx 

Jsv  Jr1 

=  ~  f  y  g'(x,y)(Cu- f)(x)dflx  + J  g'(x,y)(bu)(x)  dTx^J  (78) 

which  again  may  be  written  as  u'  =  M'  (Cu  —  /).  Note,  this  is  an  exact  formula  for  w'. 


Remarks 


1.  When  a  mesh-based  method,  such  as  finite  elements,  is  employed,  the  coarse  scales, 
u,  are  referred  to  as  the  resolved  scales ,  and  the  fine  scales,  u' ,  are  referred  to  as 
the  subgrid  scales.  The  coarse-scale  equation  is  often  referred  to  as  a  subgrid-scale 
model. 

2.  Cu  —  f  is  the  residual  of  the  resolved  scales.  It  consists  of  a  smooth  part  on  element 
interiors  (i.e. ,  f Y)  and  a  jump  term  pm]  across  element  interfaces  (i.e. ,  T'). 

3.  The  subgrid  scales  v!  are  driven  by  the  residual  of  the  resolved  scales. 

Upon  substituting  (78)  into  the  equation  for  the  coarse  scales,  (67)  is  arrived  at,  where 

(C*w,  M' (Cu  —  /))  =  -[  f  (C*w)(y)g'(x,y)(Cu- f)(x)dClxdClv 

Jn  Jo. 


[  [  (£* w)  (y) g' (x,  y)  {Cu  -  /)  (x)  dflx  dCly 
Jo.'  Jn1 

f  f  {C*w){y)g'{x,y)[bu]{x)  dTxdCly 
Jn '  Jr> 

[  I  Ww\{y)g'{x,y){Cu~  f){x)d£lxdFy 

Jv  Jn' 

I  I  (y) g'(x,  y)  [6ft]  (x)  dr x  dr y 

Jr'  Jr' 


=-  -EE 


e— i  ;= 1 


’ne  Jn1 


(C*w)(y) g'(x,  y)(Cu  -  f)(x)  dflx  dfly 


In “  Jr‘ 


I r»  Jn1 


(C*w) (y)g'(x,y)(bu)(x)  dT x  dfly 

l 

(b*w)(y) g'(x,  y){Cu  -  f)(x)  dnxdTy 


JTe  Jr1 


{b*w){y)g'(x,y){bu){x)  dT x  d£l 


(79) 
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Note,  once  again,  there  are  three  alternative  forms  due  to  the  distributional  nature  of  C*w 
and  Cu. 


Remarks 

1.  Equation  (67)  along  with  (79)  is  an  exact  equation  for  the  resolved  scales.  It  can  serve 
as  a  paradigm  for  finite  element  methods  when  unresolved  scales  are  present. 

2.  The  effect  of  the  unresolved  scales  on  the  resolved  scales  is  nonlocal. 

3.  The  necessity  of  including  jump  operator  terms  to  attain  stable  discretizations  for 
certain  problems  has  been  observed  previously  (see  Douglas  Jr.  and  Wang,  1989; 
Franca,  Hughes  and  Stenberg,  1993;  Hughes  and  Franca,  1987;  Hughes  and  Hulbert, 
1988;  Hulbert  and  Hughes,  1990;  Silvester  and  Kechkar,  1990).  The  present  result 
demonstrates  that  the  jump  operator  terms  may  be  derived  directly  from  the  governing 
equations. 

4.  Equation  (79)  illustrates  that  the  distributional  part  of  C*w  and  Cu  needs  to  be  included 
in  a  consistent  stabilized  method.  Classically,  these  terms  have  been  omitted,  which  has 
led  to  some  problems.  Jansen  et  al.  (1999)  first  observed  the  need  to  include  the  effect 
of  the  distributional  term.  In  their  approach,  rather  than  explicitly  including  the  jump 
terms,  a  variational  reconstruction  of  second-derivative  terms  is  employed.  Jansen  et  al. 
(1999)  showed  that  significant  increases  in  accuracy  are  attained  thereby.  The  method 
presented  by  Jansen  et  al.  (1999)  is  similar  to  one  presented  by  Bochev  and  Gunzburger 
(2003),  who  refer  to  procedures  of  this  kind  as  weakly  consistent. 

Equation  (67)  can  be  concisely  written  as 


B(w,  u\  g')  =  L(w;  g')  \/w  G  V 


(80) 


where 


B(w,u;g')  =  a(w,u)  +  (C*w,  M' (Cu))  (81) 

L(w-,g')  =  (w,  f)  +  (C*w,  M' f)  (82) 

Note  that  B(-,  •;  •)  is  bilinear  with  respect  to  the  first  two  arguments  and  affine  with  respect  to 
the  third  argument;  L(-;  •)  is  linear  with  respect  to  the  first  argument  and  affine  with  respect 
to  the  second.  Equations  (80)-(82)  are  valid  in  both  the  smooth  and  rough  cases,  with  the 
distributional  interpretation  appropriate  in  the  latter  case. 


Numerical  method  An  approximation ,  g'  ss  g',  is  the  key  ingredient  in  developing  a 
practical  numerical  method.  It  necessarily  entails  some  form  of  localization.  The  numerical 
method  is  written  as  follows: 


B(wh,uh ;  g')  =  L(wh ;  g')  Mw  G  V 


(83) 
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v!  and  a  posteriori  error  estimation  Note  that  u'  =  u  —  u  is  the  error  in  the  coarse 
scales.  The  formula  u'  =  M'(Cuh  —  /)  is  a  paradigm  for  a  posteriori  error  estimation.  Thus, 
it  is  plausible  that 

v!  «  M'(jCuh  -  /)  (84) 

where  M'  w  M',  is  an  a  posteriori  error  estimator,  which  can  be  used  to  estimate  coarse-scale 
error  in  any  suitable  norm,  for  example,  the  W^-norm,  0  <  p  <  oo,  0  <  s  <  oo.  (Keep  in 
mind,  Cu  is  a  Dirac  distribution.)  An  approximation  of  the  fine-scale  Green’s  function,  g'  ss  g', 
induces  an  approximation  M'  ss  M’\  see  (78).  Conversely,  an  a  posteriori  error  estimator  of 
the  form  (84)  may  be  used  to  infer  an  approximation  of  the  Green’s  function  and  to  develop 
a  numerical  method  of  the  form  (83). 

A  particularly  insightful  form  of  the  estimator  is  given  by 


Av) «  -  / 

g'(x,y)(Cuh  -f)(x)dClx  -  [  g,(a:,j/)[6'u,,](a;)dra, 

JO,’ 

Jr' 

Note  that  the  residuals  of  the  computed  coarse-scale  solution,  that  is  Cuh  —  /  and  \buh\,  are 
the  sources  of  error,  and  the  fine-scale  Green’s  function  acts  as  the  distributor  of  error. 

There  seems  to  be  agreement  in  the  literature  on  a  posteriori  estimators  that  either  one,  or 
both,  the  residuals  are  the  sources  of  error.  Where  there  seems  to  be  considerable  disagreement 
is  in  how  these  sources  are  distributed.  From  (85),  we  see  that  there  is  no  universal  solution  to 
the  question  of  what  constitutes  an  appropriate  distribution  scheme.  It  is  strongly  dependent 
on  the  particular  operator  C  through  the  fine-scale  Green’s  function.  This  result  may  serve  as 
a  context  for  understanding  differences  of  opinion  which  have  occurred  over  procedures  of  a 
posteriori  error  estimation. 


Remark 

An  advantage  of  the  variational  multiscale  method  is  that  it  comes  equipped  with  a  fine-scale 
solution  which  may  be  viewed  as  an  a  posteriori  estimate  of  the  coarse-scale  solution  error. 

Discussion 

It  is  interesting  to  examine  the  behavior  of  the  exact  counterpart  of  (85)  for  different  operators 
of  interest.  Assume  that  uh  is  piecewise  linear  in  all  cases,  and  that  /  =  0. 

First  consider  the  Laplace  operator,  C  =  —A,  b  =  d/dn.  In  this  case,  Cuh  =  0,  and  the 
interface  residual,  [6'u/1] ,  is  the  entire  source  of  error.  Keeping  in  mind  the  highly  local  nature 
of  the  Green’s  function  for  the  Laplacian,  a  local  distribution  of  would  seem  to  be  a 

reasonable  approximation.  The  same  could  be  said  for  linear  elasticity,  assuming  there  are  no 
constraints,  such  as,  for  example,  incompressibility,  or  unidirectional  inextensibility. 

Next  consider  the  advection-diffusion  operator,  C  =  a  •  V  —  kA,  =  lnduh / dnj* ,  where 

a  is  a  given  solenoidal  velocity  field,  and  k  >  0,  the  diffusivity,  is  a  positive  constant.  In  the 
case  of  diffusion  domination  (i.e.,  advective  effects  are  negligible),  C  ss  — kA,  and  the  situation 
is  the  same  as  for  the  Laplacian.  On  the  other  hand,  when  advection  dominates,  Cuh  ss  a-X/uh, 


*This  follows  from  the  continuity  of  advective  flux. 
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and  \buh ]  =  \ndvr /dri\  may  be  ignored.  This  time,  the  element  residual,  £uh,  is  the  primary 
source  of  error.  A  local  distribution  scheme  would  seem  less  than  optimal  because  the  Green’s 
function  propagates  information  along  the  integral  curves  of  —a  (keep  in  mind  that  the  Green’s 
function  is  for  the  adjoint  operator,  £*),  with  little  amplitude  decay.  This  means  that  there  is 
an  approximately  constant  trajectory  of  error  corresponding  to  the  residual  error  £uh,  in  the 
element  in  question. 

Finally,  consider  the  Helmholtz  operator,  £  =  —A  —  fc2,  b  =  d/dn,  where  fc  is  the  wave 
number.  If  k  is  real,  we  have  propagating  waves ,  whereas  if  k  is  imaginary,  we  have 
evanescent  (decaying)  waves.  In  the  latter  case,  the  Green’s  function  is  highly  localized; 
as  fc  — »■  0  the  Green’s  function  approaches  that  for  the  Laplacian,  as  |fc|  — >  oo  the  Green’s 
function  approaches  —k~26,  a  delta  function.  In  the  case  of  propagating  waves,  the  Green’s 
function  is  oscillatory.  In  general,  for  |fc|  large,  the  dominant  source  of  error  is  the  element 
residual,  £uh  =  — k2uh .  As  |fc|  — >  0,  the  interface  residual,  pra71]  =  \duh/dn\ ,  dominates. 

3.3.  Hierarchical  p-refinement  and  bubbles 

Hierarchical  p-refinement  plays  an  important  role  in  clarifying  the  nature  of  the  fine-scale 
Green’s  function,  g',  and  provides  a  framework  for  its  approximation.  Some  notations  are 
required.  Let 

^nodes 

uh  =  Naua  (likewise  wh)  (86) 

A—l 

where  Na  is  a  finite  element  shape  function  associated  with  the  primary  nodes,  A  = 
1, 2,  ■  ■  ■  ,  hnodes,  and  ua  is  the  corresponding  nodal  value;  and  let 

^nodes 

u  =  ^  N'au'a  (likewise  w')  (87) 

A-l 

where  N'A  is  a  hierarchical  finite  element  shape  function  associated  with  the  additional  nodes, 

A  =  1,2,- ••  ,n'nodes,  and  u'a  are  the  corresponding  hierarchical  degrees  of  freedom.  For 

example,  let  uh  be  expanded  in  piecewise  linear  basis  functions  and  v!  in  hierarchical  cubics 
(see  Fig.  20).  Note,  bubble  functions  are  zero  on  element  boundaries.  An  illustration  in  one 
dimension  is  presented  in  Figure  21. 

Substituting  (86)  and  (87)  into  (57)-(60),  and  eliminating  u'A  by  static  condensation 
results  in 

B(wh,uh ;  g')  =  L(wh;  g')  \/wh  G  V  (88) 

where 

B(wh,  uh\  g')  =  a{wh,uh)  +  (C*wh,M\Cuh))  (89) 

L(wh;  §')  =  (wh,f)  +  (C*wh,M'f)  (90) 

and 

(C*wh,M\Cuh))  =  -  [  [  (£*wh)(y)g\x,y)(£uh)(x)dnxday  (91) 

JQ  JQ 
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uh  =  standard  linears  (•) 
u'  =  hierarchical  cubics  (o) 


Figure  20.  Hierarchical  cubics  in  two  dimensions. 
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standard  linear 
shape  functions 


quadratic  bubbles 


cubic  bubbles 


etc. 


Figure  21.  Finite  element  shape  functions  and  polynomial  “bubbles”. 


where 


g '(*,»)=  E  kmUk")-1}  n'b(x) 

J  Ad 


A,B— 1 


k"  =  [K"AB\ 

KAb  =  a(NA,N'B) 


(92) 


(93) 

(94) 


Remarks 

1.  Recall,  Cuh  and  C*wh  are  Dirac  distributions  in  the  finite  element  case  (cf.  (91)  and 
(79)). 

2.  Hierarchical  p-refinement  generates  an  approximate  fine-scale  Green’s  function,  g'  « 
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3.  For  implementational  purposes ,  it  is  more  convenient  to  rewrite  (89)-(91)  in  forms 
avoiding  Dirac  distributions.  This  can  be  accomplished  by  using  the  integration- by-parts 
formulas,  viz., 

^nodes 

(C*wh,  M'(Cuh  —  /))  =  —  £  a(wh,N'A)  (. K ")~l  (a(N'B,uh)  -  (N’B,  /))  (95) 

z — '  L  J  AB 

A,B=1 

which  amounts  to  the  usual  static  condensation  algorithm. 

4.  A  posteriori  error  estimation  for  the  coarse-scale  solution,  uh,  is  provided  by  the  fine- 
scale  solution  (see  (78)  and  (84)): 


^nodes 

Av)  =  ~  E  NA(y)[(K")-1]AB{a(NlB,uh)-(N'B,f))  (97) 

A,B= 1 


The  quality  of  this  estimator  depends  on  the  ability  of  {NlA}nACAles  to  approximate  the 
fine  scales,  or  equivalently,  the  quality  of  the  approximation  g'  ss  g' . 

5.  Note  that  the  fine-scale  Green’s  function  only  depends  on  the  hierarchical  basis  (see 
(92)-(94)).  The  exact  fine-scale  Green’s  function  corresponds  to  the  limit  p  — >  oo. 

6.  The  fine-scale  Green’s  function  is  nonlocal ,  but  it  is  computed  from  a  basis  of  functions 
having  compact  support.  For  example,  in  the  two-dimensional  case,  the  basis  consists 
of  bubbles,  supported  by  individual  elements,  and  edge  functions,  supported  by  pairs 
of  elements  sharing  an  edge.  The  three-dimensional  case  is  similar,  but  somewhat  more 
complicated;  the  basis  consists  of  bubbles,  face  and  edge  functions.  In  three  dimensions, 
pairs  of  elements  support  face  functions  whereas  the  number  of  elements  supporting 
edge  functions  depends  on  the  topology  of  the  mesh. 

7.  In  two  dimensions,  by  virtue  of  the  convergence  of  hierarchical  p-refinement,  the  exact 
fine-scale  solution  may  be  decomposed  into  a  finite  number  of  limit  functions  -  one 
bubble  for  each  element  and  one  edge  function  for  each  pair  of  elements  sharing  an 
edge.  In  one  dimension  the  situation  is  simpler  in  that  only  bubbles  are  required.  The 
three-dimensional  case  is  more  complex  in  that  bubbles,  face  and  edge  functions  are 
required. 
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8.  Polynomial  bubbles  are  typically  ineffective,  but  so-called  residual-free  bubbles 
(Brezzi  et  al.  (1997))  are  equivalent  to  exactly  calculating  element  Green’s  functions. 
This  approximation  works  exceptionally  well  in  some  important  cases  and  will  be 
discussed  in  more  detail  later  on. 

9.  If  only  bubble  functions  are  considered  in  the  refinement,  coupling  between  different 
elements  and  all  jump  terms  is  eliminated.  In  this  case, 


(. C.*wh,M'(Cuh ))  =  - 


(£*wh)(y) g'(x,  y)(£uh)(x )  dflx  dfly 


(98) 


Cl'  J  Cl' 


where  17'  has  replaced  f 7,  and  now  the  effect  of  u'  is  nonlocal  only  within  each  element. 
The  approximate  Green’s  function  g',  is  defined  element-wise  and  takes  on  zero  values 
on  element  boundaries. 

10.  It  may  be  observed  that  the  fine-Green’s  function  formula  has  an  interpretation 

analogous  to  projection  methods  in  linear  algebraic  systems,  such  as,  for  example, 
multigrid  methods.  To  explicate  this  analogy,  the  notation  of  Saad  (1995),  Chapter  5,  is 
adopted.  Let  1C  be  an  m-dimensional  subspace  of  ffi" ,  and  let  A  be  an  n  x  n  real  matrix. 
The  objective  is  to  solve  Ax  =  b  for  x  £  M” ,  where  b  £  Wl  is  a  given  vector.  Assume  an 
initial  guess,  xq  £  ffi” ,  and  determine  an  approximate  solution  x  €  Xq  +  /C,  such  that  the 
residual  f  =  b—Ax  _L  1C.  If  x  is  written  as  x  =  Xo+<5,  where  6  £  1C,  and  ro  =  b—Ax o,  then 
b  —  A(x o  +  S)  _L  1C,  or  equivalently,  ro  —  AS  _L  1C.  In  terms  of  a  basis  {vi,V2,  •  ■  • ,  vm} 
of  /C,  the  approximate  solution  becomes  x  =  Xq  +  Vy,  where  V  =  [iq,  V2,  ■  ■  ■ ,  vm\, 
an  n  x  m  matrix.  The  orthogonality  condition  becomes  VT AVy  =  VTrg,  and  thus 
S  =  x  —  x o  =  P(PTAP)_1PTro,  assuming  nonsingularity  of  VT  AV.  Note  the  similarity 
of  this  result  to  that  of  (92)  and  (97):  6  ~  u' ,  x  ~  u,  xq  ~  uh ,  V  ~  [N[,N2,  ■  ■  ■ ,  N',  ], 

VTAV  ~  K”,  VTr0  ~  [a(IV{,#)  -  (iV{ ,/),..., a{N'n,  ,uh)  -  {N'n,  ,/)]T,  and 

nodes  nodes 

V{VT AV)~1VT  ~  g' .  In  multilevel  solution  strategies,  the  fine-scale  space  here  may 
be  analogized  to  the  “coarse  grid”,  and  S  is  obtained  by  restriction  to  the  coarse-grid 
subspace  (i.e.,  VTr0),  a  coarse-grid  correction  (i.e.,  VT AVy  =  VTr0),  and  prolongation, 
or  interpolation  (i.e.,  Vy). 

11.  Defect- correction  techniques  have  become  popular  in  the  multigrid  community  to 
obtain  stable  second-order  accurate  solutions  of  the  convection-diffusion  equation. 
These  methods  use  upwind  schemes  for  relaxation  (stability)  and  central  difference 
schemes  for  residual  evaluation  (accuracy).  A  standard  reference  is  Hemker  (1981); 
see  also  Trottenberg,  Oosterlee  and  Schuller  (2001).  More  generally,  defect-correction 
is  a  powerful  abstraction  for  various  iterative  methods,  including  Newton,  multigrid, 
and  domain  decomposition.  It  can  also  be  extended  to  multiscale  problems.  See  Lai 
(1981)  for  an  example  of  how  defect-correction  is  used  in  the  treatment  of  small-scale 
fluctuations  in  aeroacoustics. 


3-4-  Residual-free  bubbles 

The  concept  of  residual-free  bubbles  has  been  developed  and  explored  in  Baiocchi,  Brezzi 
and  Franca  (1993),  Brezzi  et  al.  (1997),  Brezzi  and  Russo  (1994),  Franca  and  Russo  (1996), 
Russo  (1996a,  1996b).  The  basic  idea  is  to  solve  the  fine-scale  equation  on  individual  elements 
with  zero  Dirichlet  boundary  conditions.  For  example,  the  objective  is  to  find  u'  £  V',  such 
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that  Vw£  V, 


{II')*  Cu'  =  -(IT )*(Cu-f)  on  fle 

u'  =  0  on  re 


e  =  1, 2,  ■ 


>  Tie  1 


(99) 


Noting  that  u  can  be  expressed  in  terms  of  the  coarse-scale  basis  having  support  in  the  element 
in  question,  a  fine-scale  basis  of  residual-free  bubbles  can  be  constructed  for  each  element,  i.e. 


(lT)*CN'a  =  -(n’)*(£Na  -  f)  on  IT  \ 

N'a  =  0  onTe  J 


e=  1,2,  -  ,nei  (100) 


where  a  =  1,2,  ■■■  ,nen  is  the  local  numbering  of  the  primary  nodes  of  element  e.  Thus,  to 
each  coarse-scale  basis  function  Na,  solve  (100)  for  a  corresponding  residual-free  bubble  N’a. 
Consequently,  the  maximal  dimension  of  the  space  of  residual-free  bubbles  for  element  e  is 
nen.  It  is  typical,  however,  that  the  dimension  is  less  than  nen. 

Brezzi  and  Russo  (1994)  have  constructed  residual-free  bubbles  for  the  homogeneous 
advection-diffusion  equation  assuming  the  coarse-scale  basis  consists  of  continuous,  piecewise 
linears  on  triangles.  For  this  case  nen  =  3,  but  the  dimension  of  the  space  of  residual-free 
bubbles  is  only  one.  Let  Be  denote  the  residual- free  bubble  basis  solution  of  the  following 
problem: 


CBe  =  1  on  fie  1 

Be  =  0  onTe  / 


(101) 


Note  that,  due  to  the  fact  the  coarse-scale  space  consist  only  of  piecewise  linears,  combined  with 
the  fact  that  the  fine-scale  space  satisfies  zero  Dirichlet  boundary  conditions,  the  projection 
operator,  IT,  present  in  the  general  case,  namely  (100),  can  be  omitted  and  (101)  can  be  solved 
in  the  strong  sense.  However,  in  order  to  avoid  potential  linear  dependencies,  in  general,  (100) 
needs  to  be  respected. 

In  principle,  the  computation  of  the  residual-free  bubble  should  involve  the  solution  of 
one  or  more  partial  differential  equations  in  each  element.  This,  however,  can  be  done  in  an 
approximate  way,  using  a  suitable  subgrid  in  each  element,  as  in  Brezzi,  Marini  and  Russo 
(1998)  and  Franca,  Nesliturk  and  Stynes  (1998).  This  leads  to  the  idea  of  the  stabilizing 
subgrid:  if  one  does  not  eliminate  the  discrete  bubbles,  one  can  think  of  solving  the  Galerkin 
method  on  the  enriched  subgrid,  having  the  beneficial  effects  of  the  bubble  stabilization.  See 
Brezzi  and  Marini  (2002),  Brezzi  et  al.  (2003),  and  Brezzi,  Marini  and  Russo  (2004)  for  further 
development  of  this  idea. 


3.5.  Element  Green’s  functions 


The  idea  of  employing  an  element  Green’s  function  was  proposed  in  the  initial  work  on 
the  variational  multiscale  method  (Hughes  (1995)).  In  place  of  (63)-(64),  the  Green’s  function 
problem  for  each  element  is  solved: 


(H')*c*g'e(x,y)  =  \/x  €  1 

g'e(x,y)  =  0  Vxde  J 


e=  1,2,  •  ,nel  (102) 


Use  of  element  Green’s  functions  in  place  of  the  global  Green’s  function  amounts  to  a  local 
approximation , 


g'(x,y)  =  g'e(x,y)  Vx,  y  €  fie,  e  =  1, 2,  •  •  •  ,  nei 


(103) 
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The  upshot  is  that  the  subgrid  scales  vanish  on  element  boundaries,  i.e. , 

v!  =  0  on  Te,  e  =  1,  2,  ■  ■  ■  ,nei 


(104) 


This  means  the  subgrid  scales  are  completely  confined  within  element  interiors. 

There  is  an  intimate  link  between  element  Green’s  functions  and  residual-free  bubbles.  This 
idea  was  first  explored  in  Brezzi  et  al.  (1997),  in  which  it  was  shown  that,  for  the  case  governed 
by  (101), 


Be{y)=  [  g 'e{x,y)dttx 
J  ne 


(105) 


This  result  can  be  easily  derived  as  follows: 

[  g 'e(x,y)dClx  =  (g',l)f 

J 


=  (g'e,£-Be)nS 
=  «(g  e,Beh‘x 
=  (C*g'e,Be)ne 
=  (S,Be)ae 
=  Be(y) 


(106) 


Another  way  to  derive  (105)  is  to  appeal  to  the  general  formula  for  the  Green’s  function  in 
terms  of  a  fine-scale  basis,  namely  (92).  Specialized  to  the  present  case,  (92)  becomes 


Note  that 


g'e(x,  y )  =  Be(x)  (a(Be,  Be)ne)  1  Be(y) 

a(Be:Be)ne  =  (Be,CBe)ne 
=  (Be,l)n* 

=  [  Bed.n 

The  result  follows  by  integrating  (107)  and  using  (108), 

[  g'e(x,y)dnx  =  f  Be{x)dQ.x{a(yBe,Be)o.e)~1  Be{y) 
Jne  Jne 

=  Be(y) 


(107) 


(108) 


(109) 


Remarks 


1.  In  general,  the  relationship  between  a  fine-scale  basis  and  a  Green’s  function  is  given 
by  (92).  The  result  (105)  is  special  to  a  residual-free  bubble  governed  by  (101). 

2.  In  general,  the  coarse-scale  residual  will  not  be  constant  within  an  element.  For  example, 
suppose  the  residual  is  a  linear  polynomial  in  x.  Then  there  are  two  residual-free  bubbles, 
one  corresponding  to  1  and  one  corresponding  to  x,  say  and  B^\  respectively.  B^ 
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is  the  same  as  Be,  and  is  given  by  (105).  may  be  determined  in  the  same  way  as 
(106): 


g'e(x,y)xdflx 


(g'e.^Oj 

(CCB^)ne 

a(g'e,B^)n% 

(C*g'e,BW)n% 

(S,B^)ae 

Bi'Hy) 


(110) 


As  may  be  seen,  B ^  is  the  first  moment  of  g'e.  This  is  the  general  case:  If  the  coarse- 
scale  residual  is  expanded  in  terms  of  a  basis  of  functions,  then  the  residual-free  bubbles 
are  the  moments  of  g'e  with  respect  to  the  elements  of  that  basis. 


3.6.  Stabilized  methods 


Classical  stabilized  methods  are  generalized  Galerkin  methods  of  the  form 


a(wh,  uh)  +  (L wh,  r(Cuh  -  /)) n'  =  ( wh ,  /) 
where  L  is  typically  a  differential  operator ,  such  as 

(111) 

+ 

II 

J 

Galer kin/least-squares  (GLS) 

(112) 

L  —  +-^adv 

SUPG 

(113) 

* 

1 

II 

J 

Multiscale 

(114) 

and  r  is  typically  an  algebraic  operator.  SUPG  is  a  method  defined  for  advective-diffusive 
operators,  that  is,  ones  decomposable  into  advective  (£adv)  and  diffusive  (jCdiff)  parts.  A 
stabilized  method  of  the  form  (114)  is  referred  to  as  a  “multiscale”  stabilized  method  for 
reasons  that  will  become  apparent  shortly. 


3.6.1.  Relationship  of  stabilized  methods  with  subgrid-scale  models  It  was  shown  in  Hughes 
(1995)  that  a  stabilized  method  of  multiscale  type  is  an  approximate  subgrid-scale  model 
in  which  the  algebraic  operator  r  approximates  the  integral  operator  M'  based  on  element 
Green’s  functions , 

r  =  —M'  «  —  M'  (115) 


Equivalently, 


T-S(y-x)=  g'e(x,  y)  k.  g'(x,  y) 


(116) 


The  result  follows  from  the  calculation 


[  [  (-C*wh)(y)g'e(x,y)(Cuh  -  f)(x)dflxdfly 

J  O'  -hr 

=  [  [  {-C*wh){y)T  ■  S(y  -  x){Cuh  -  f)(x)  dflx  dCl% 

■hi’  .hi' 

=  [  (- C*wh)(x)r  ■  ( Cuh  -  f)(x)  dtlx 

Js V 


(117) 
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3.6.2.  Formula  for  r  based  on  the  element  Green’s  function  The  approximation 

T-S{y-x )  «g 'e(x,y) 

suggests  defining  r  by 

/  /  r  •  <5(y  —  a;)  =  /  /  gg(:r,  it)  cttda;  cftlj, 

Jne  Jne  Jne  Jne 


meas(fie)  Jn. 


y)  dFlx  d^ly 


Remarks 


(118) 

(119) 

(120) 


1.  The  element  mean  value  of  the  Green’s  function  provides  the  simplest  definition  of  r. 

2.  This  formula  is  adequate  for  low-order  methods  (h-adaptivity).  For  higher-order 
methods  (p-adaptivity),  accounting  for  variation  of  r  over  an  element  may  be  required. 
In  this  case,  it  may  be  assumed,  for  example,  that  r  =  t{x ,  y)  is  a  polynomial  of 
sufficiently  high  degree.  Given  an  element  Green’s  function  g’e,  an  equivalent  function  r 
can,  in  principle,  always  be  calculated.  Consequently,  there  is  a  generalized  stabilized 
method ,  that  is,  a  method  of  the  form, 

a(wh,uh)  —  ’S~'  [  [  C*wh(y)T(x,y)(Cuh  -  f)(x)  dnxdQy  =  (wh,  f)  (121) 

e  Jne  J 

equivalent  to  the  element  Green’s  function  method.  The  generalized  stabilized  method 
involves  determining  r  such  that  the  following  equivalence  condition  is  satisfied 


C*wh(y)T{x,y)(Cuh  -  f)(x)  dVLxdFLy 
C*wh(y)g'e(x,y)(Cuh  -  f)(x)  dflx  dFlv 


\/wh,uh  €  V  (122) 


Jne  J Qe 

Thus  a  full  equivalence  exists  as  indicated  in  Figure  22.  Examples  of  the  calculation  of 
r  by  way  of  this  procedure  will  be  presented  subsequently. 


Examples 

Two  one-dimensional  examples  are  considered,  the  advection-diffusion  equation  and  the 
Helmholtz  equation.  In  both  cases  it  is  assumed  that  standard  piecewise  linears  are  used  for 
the  coarse-scale  basis.  Consequently,  determination  of  the  element  Green’s  function  may  be 
performed  using  the  strong-form  counterpart  of  (102),  namely 


£*g  10m/)  =  S(x-y)  Vxefle  ) 

g'e(x,y)  =  0  VxeTe  J 


e  —  1 , 2 ,  *  ,  ne  i 


(123) 


Remark 

In  one  dimension ,  use  of  the  element  Green’s  function  g’e  results  in  u  =  u  +  v!  being 
pointwise  exact  for  any 

Anodes 

U=  Y,  Naua  (124) 

A—l 

and  u  being  an  end-node  interpolant  of  u.  This  result  holds  for  all  problems. 
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residual  free  bubbles  <©  ©>  element  Green’s  functions 


stabilized  methods 


Figure  22.  Stabilized  methods  can  be  constructed  which  are  equivalent  to  methods  based  on  element 
Green’s  functions,  which  in  turn  are  equivalent  to  methods  developed  from  the  residual-free  bubbles 
concept.  Historically,  stabilized  methods  preceded  the  residual-free  bubbles  and  element  Green’s 
function  approaches,  and  there  are  certain  stabilized  methods  (e.g.,  SUPG  and  GLS)  which  have  been 
established  for  some  time  that  are  not,  strictly  speaking,  equivalent  to  these  concepts.  Nevertheless, 
they  have  been  justified  independently  by  mathematical  analysis  and  numerical  testing,  and,  in  the 
case  of  SUPG,  may  be  viewed  as  a  simplified  variational  multiscale  method. 


Figure  23.  Element  Green’s  function  for  the  one-dimensional  advection-diffusion  equation  as  a  function 

of  element  Peclet  number  (a). 
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oo 


a  >  1 

advection  dominated 


a  <  1 

diffusion  dominated 


a 


0 


Figure  24.  Schematic  behavior  of  the  element  Green’s  function  for  the  advection-diffusion  operator, 

for  y  fixed. 


Advection-diffusion  equation 

Let 

—  -^adv  H"~  -^diff 

(125) 

where 

-^adv  —  tl  ~~r 
ax 

(126) 

r  d 2 

£diff  “  ^ dx 2 

(127) 

and  a  and  k  are  assumed  to  be  positive  constants.  Consider  the  homogeneous  Dirichlet 
problem 


Cu  =  f 

in  CL  =  [0,  L\ 

(128) 

u  =  0 

on  T  =  {0,  L} 

(129) 

C*  = 

r*  _i_  r* 

•^adv  '  ^-'diff 

(130) 

= 

^adv  -^difF 

(131) 

The  solution  of  (123)  is  given  by 


g'eO v,y) 


Ci(y)(  l-e-2“f) 
C2(y){e~2a %  -  e~2a) 


x  <  y 
x  >  y 


(132) 
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where 


Ci(y)  = 

l-e-M1"*) 

(133) 

a  (1  —  e~2a) 
e2<  -  1 

C2(y)  = 

(134) 

a  (1  —  e~2a) 

a  = 

—  (element  Peclet  number) 

2k 

(135) 

and  h  =  meas(fle)  is  the  element  length.  This  element  Green’s  function  is  shown  in  Figure  23 
for  various  element  Peclet  numbers.  For  fixed  y  it  behaves  as  shown  in  Figure  24. 

By  virtue  of  the  fact  that  the  coarse  scales  are  piecewise  linear  and  Cuh  and  C*wh  are 
therefore  constants  on  each  element,  the  simple  formula  (120)  suffices  to  define  a  constant  r 
satisfying  (122),  that  is,  one  equivalent  to  use  of  the  element  Green’s  function,  viz., 


T 


h  (  ,  1\ 

2a  \  a  J 


(136) 


Remarks 


1.  This  is  a  well  known  formula  from  the  theory  of  stabilized  methods.  It  was  originally 
derived  using  Fourier  methods  on  regular  meshes  and  assuming  constant  coefficients. 

2.  This  t  results  in  a  nodally  exact  stabilized  method  for  piecewise  linear  Na ’s  and 
element-wise  constant  a,  k  and  /.  Remarkably,  element  lengths  need  not  be  uniform. 

3.  By  virtue  of  (105)  and  (135),  the  same  optimal  r  is  obtained  from  the  residual- free 
bubble,  namely, 


T 


1 

meas(fle) 

h  ,  , 

—  (coth  a 

2  a 


Be(x )  dii 


J 

-h 


(137) 


4.  If  the  fine-scales  are  modeled  by  a  single  quadratic  polynomial  bubble  in  each  element, 
the  approximate  element  Green’s  function  is  given  by 


g  e(x>y) 


N'(x)N'(y) 

a(N',N') 


(138) 


Employing  (120)  results  in 


h 2  h  a 

12k  2 a  3 


(139) 


Comparing  (139)  with  (137)  reveals  that  t  is  much  larger  than  r  in  advection  dominated 
cases.  See  Figure  25.  This  results  in  a  method  for  the  coarse  scales  which  is  overly 
diffusive. 

Although  the  addition  of  the  quadratic  polynomial  bubble  produces  an  overly  diffuse 
method,  it  is  a  curious  fact  that  in  conjunction  with  the  addition  of  an  appropriately 
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Figure  25.  Behavior  t  and  t  as  a-*  oo. 


defined  artificial  diffusivity  in  the  fine-scale  equation,  the  modified  method  is  capable 
of  generating  the  optimal  value  of  r  given  by  (137).  This  may  be  seen  as  follows:  in  the 
fine-scale  equation,  (59),  replace  the  term  a(wfu')  by  a(w' ,u')  +  a'(w',ii')  where 

a'(w',  u,')  =  (Vw',  k'Vu')  (140) 


in  which  u'  is  the  fine-scale  artificial  diffusivity.  In  the  present,  one-dimensional 
case,  this  term  is  simply  K'^t)  and  reduces  to  element  integrals  of  the  derivatives 
of  the  quadratic  polynomial  bubbles.  From  the  previous  development  it  is  clear  that 
this  will  produce  a  r  of  the  same  form  as  (139)  but  with  k  replaced  by  k  +  k',  namely 


r 


7 


h2 

12(k+k') 


(141) 


Equating  t'  with  r  leads  to  the  following  definition  for  k': 

The  “effective  diffusivity”  in  the  fine  scales  is 


ah 


k  +  u‘  =  —  =  Oih) 

6 


(142) 


(143) 


This  method  is  identical  to  the  one  produced  by  the  exact  element  Green’s  function  and 
thus  is  superconvergent  with  respect  to  the  nodal  values  (i.e. ,  it  is  exact  at  the  nodes), 
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and  attains  optimal  coarse-scale  convergence  rates  in  the  L2  and  H1  norms  (i.e.,  0(h 2) 
and  O(h),  resp.).  It  is  a  somewhat  remarkable  fact  that  adding  artificial  diffusivity  to  the 
fine-scale  equation  actually  reduces  the  diffusive  effect  in  the  coarse-scale  equation.  If 
an  0(h)  effective  diffusivity  of  the  form  (143)  was  added  to  all  scales  (i.e.,  to  (56)),  then 
the  method  would  amount  to  the  classical  artificial  diffusivity  method  which  is  overly 
diffusive  and  limited  to  0(h)  in  L2  and  0(h1^2)  in  H1,  independent  of  the  polynomial 
order  of  the  elements  employed.  It  may  be  concluded  that  fine-scale  artificial  diffusivity 
is  a  viable  modeling  tool  in  multiscale  analysis  and  one  that  is  very  different  than, 
and  superior  to,  classical  artificial  diffusivity  in  all  scales.  Of  course,  to  obtain  (142), 
the  optimal  value  of  r,  given  by  (137),  needed  to  be  known.  (The  idea  of  a  fine-scale 
artificial  diffusivity  is  due  to  Guermond  (2001),  and  has  been  further  analysed  in  Brezzi 
et  al.  (2000).) 

5.  The  approximate  fine-scale  solution  is  given  by 

u'(y)  «  -  [  g'(x,y)(Cuh  -f)(x)dflx  -  [  g' (x,y)\buh\(x)  dT x  (144) 

h v  Jt' 

where  g'  is  an  approximation  to  the  fine-scale  Green’s  function,  g' .  For  the 
multidimensional  advection-dominated  case,  the  jump  term  can  be  neglected  and  the 
approximation  Cuh  «  a  ■  X/uh  can  be  used.  Furthermore,  employing  the  approximation 
(116)  results  in 

u'  «  — r(a  •  Vuh  —  f)  (145) 

Assuming  r  has  the  form  h/(2\a\)  in  the  advection-dominated  case,  (145)  becomes 

M'W~2R(a' WW)  (146) 

This  may  be  used  as  a  simple,  local,  error  estimator  for  the  advection-dominated  case. 


Helmholtz  equation  The  set-up  is  identical  to  the  previous  example.  Consider  the 
Dirichlet  problem,  (128)-(129),  for  the  Helmholtz  equation  in  which 


C  = 


dx 2 


k 2 


(147) 


is  the  Helmholtz  operator  and  k  is  the  wave  number.  Assume  k  G  R,  corresponding  to  the 
case  of  propagating  waves.  Note  C*  =  C.  The  Green’s  function  for  an  element  of  length  h 
is  given  by 


g'e(tOl) 


sm(kh(l  +  $) / 2)  sin(fch(l  —  rj) /2) 

ksin(kh)  1  <  ^ 

sin(fc/i(l  —  £) /2)  sm(kh(l  +  g) /2) 

ksm(kh)  ’  ^  1 


(148) 


where  £  and  g  are  normalized,  bi-unit  coordinates, 


-1  <  id  <  +1 


(149) 


The  element  Green’s  function  is  depicted  in  Figure  26  for  various  element  wave  numbers. 
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(c)  kh  =  2tx  (d)  kh  =  4-rr 

Figure  26.  Element  Green’s  function  ( kg'e )  for  the  one-dimensional  Helmholtz  equation  as  a  function 

of  wave  number  (kh). 


This  time  Cuh  and  C*wh  vary  linearly  over  each  element  and  satisfaction  of  the  equivalence 
condition,  (122),  entails  a  non-constant  r,  as  follows: 

t(£,v)  =T00  +  Tll^«g(,(^??)  (150) 


where 


£Vg  e(^v)d£dri 


l-l  ,7  —  1 

/*+!  /*+! 


F'rpdtdn 


Remarks 


1.  Oberai  and  Pinsky  (1998)  have  also  derived  an  element  Green’s  function  in  two 
dimensions  for  a  bilinear  rectangle  and  an  equivalent  r.  Similar  results  for  three 
dimensions  are  also  easily  derived. 
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2.  As  in  the  previous  example,  this  r  results  in  a  nodally  exact  stabilized  method  for 
piecewise  linear  Na’s  and  element-wise  constant  k  and  /,  for  an  arbitrary  non-uniform 
mesh. 

3. 7.  Summary 

In  this  section,  the  variational  multiscale  formulation  for  an  abstract  Dirichlet  problem  was 
described.  The  case  in  which  finite  elements  are  used  to  represent  the  coarse  scales  was 
considered  and  an  exact  equation  governing  the  coarse  scales  was  derived.  This  equation 
amounts  to  the  Galerkin  formulation  plus  an  additional  term  which  depends  on  the 
distributional  form  of  the  residual  (i.e.,  element  interior  and  interface  jump  terms)  and  a 
fine-scale  Green’s  function.  A  representation  was  also  derived  for  the  fine-scale  solution,  which 
amounts  to  the  error  in  the  coarse  scales.  The  results  serve  as  paradigms  for  subgrid-scale 
models  and  a  posteriori  error  estimators. 

To  understand  the  nature  of  the  fine-scale  Green’s  function,  hierarchical  p-refinement  and 
bubbles  were  considered.  An  explicit  formula  for  the  fine-scale  Green’s  function  was  obtained 
in  terms  of  the  hierarchical  (i.e.,  fine-scale)  basis.  This  formula  suggests  that  the  fine-scale 
Green’s  function  can  be  represented  in  terms  of  a  basis  of  functions  having  local  support. 
Subsequent  discussion  dealt  with  developing  practical  approximations. 

The  concepts  of  residual-free  bubbles,  element  Green’s  functions,  and  stabilized  methods 
were  reviewed,  and  relationships  among  them  were  summarized. 


4.  Space-time  Formulations 


In  Sections  2  and  3,  a  methodology  was  described  to  address  problems  in  which  scales  are 
present  that  are  viewed  as  numerically  unresolvable,  and  accurate  computation  of  the  resolvable 
scales  necessitates  incorporating  the  effects  of  the  unresolvable  scales.  The  analysis  of  Section 
3  was  restricted  to  steady  phenomena.  In  this  section  the  study  of  the  variational  multiscale 
method  is  continued,  generalizing  the  development  of  the  subgrid-scale  models  to  the  time- 
dependent  case.  A  first-order  in  time,  second-order  in  space,  non-symmetric,  linear  partial 
differential  equation  is  considered.  Creation  of  the  subgrid-scale  model  of  the  initial/boundary- 
value  problem  requires  solution  of  an  element  problem  for  the  unresolved  scales  which  is 
overspecified.  This  theoretical  impediment  is  overcome  by  way  of  an  elliptic  regularization 
procedure,  giving  rise  to  an  element  Green’s  function  problem.  In  the  limiting  case  of  the 
regularization  parameter,  the  desired  causal  Green’s  function  of  the  adjoint  operator  is  derived. 
This  enables  the  subgrid  scales  to  be  determined  analytically  and  represented  in  terms  of  the 
residual  of  the  resolved  scales. 

The  variational  setting  for  the  developments  is  space-time.  That  is,  a  fully-discrete  method 
is  developed  entirely  from  finite  element  concepts.  It  is  believed  this  provides  the  most 
theoretically  coherent  framework  for  the  derivations. 

As  in  Section  3,  the  subgrid-scale  model,  when  numerically  approximated  by  Galerkin’s 
method,  provides,  on  the  one  hand,  a  paradigm  for  bubble  function  finite  element  methods, 
and,  on  the  other  hand,  an  expose  of  the  theoretical  foundations  of  stabilized  methods.  Both 
bubbles  and  stabilized  methods  are  thus  identified  as  approximate  subgrid-scale  models.  In 
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Figure  27.  Finite  element  discretization  in  space-time. 


the  case  of  stabilized  methods,  this  identification  leads  to  formulas  enabling  the  calculation  of 
r,  the  parameter  present  in  stabilized  methods. 

To  illustrate  the  use  of  the  theory  in  calculating  r,  the  simple  case  of  an  ordinary  differential 
equation  in  time  is  considered.  The  known  optimal  value  is  arrived  at,  complementing  the 
analogous  calculation  for  the  advection-diffusion  operator  in  space,  presented  in  Section  3. 

4-1-  Finite  elements  in  space-time 

Consider  a  discretization  of  the  space-time  domain  in  question  into  time  slabs  which  are  in 
turn  discretized  into  space-time  elements.  A  schematic  illustration  is  presented  in  Figure  27. 
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Slab  Qn  is  the  space-time  domain  bounded  by  spatial  hypersurfaces  at  times  tn  and  tn+\. 
The  lateral  boundary  of  Qn  is  denoted  by  Pn.  A  decomposition  of  this  type  is  a  suitable 
starting  point  for  the  variational  formulation  of  equations  of  evolution.  A  discrete  problem 
is  solved  on  each  slab  sequentially,  starting  with  the  first.  This  gives  rise  to  a  fully-discrete 
time-stepping  algorithm.  For  appropriately  formulated  finite  element  methods,  the  space-time 
approach  is  amenable  to  a  priori  and  a  posteriori  error  estimation,  adaptive  strategies,  and 
relatively  simple  moving  domain  procedures  (see,  e.g.,  Aliabadi  and  Tezduyar  (1993),  Hughes, 
Franca  and  Hulbert  (1989),  Hughes  and  Hulbert  (1988),  Johnson  (1986,1987,1992),  Johnson 
and  Hansbo  (1992),  Johnson,  Navert  and  Pitkaranta  (1984),  Mittal  and  Tezduyar  (1994), 
Tezduyar,  Behr  and  Liou  (1992),  Tezduyar  et  al.  (1992)). 

4-2.  Subgrid-scale  modeling 

To  simplify  subsequent  writing,  let  Q  denote  a  generic  space-time  slab  bounded  at  t  =  0  by 
Oo  and  at  t  =  T  by  f It-  Thus  Q  may  represent  any  slab  or  the  entire  space-time  domain.  The 
following  definitions  are  important  for  subsequent  developments: 


Q'  = 

riel 

U«” 

e— 1 

(element  interiors) 

(152) 

p'  = 

ne  l 

l Jr* 

— 1 

(element  boundaries) 

(153) 

Q  = 

Q  =  closure(<5') 

(154) 

p  = 

dQ\(n0 

u 

(155) 

ne i  is  the  number  of  elements  in  slab  Q.  Consider  an  overlapping  sum  decomposition  in  which 
u  represents  the  resolvable  scales  and  v!  represents  the  unresolvable ,  or  subgrid ,  scales. 

^.3.  Initial/boundary-value  problem 

Consider  the  following  initial/boundary- value  problem:  Given  /  :  Q  — >  ffi,  g  :  P  — >  ffi,  and 
it(0-)  :  fio  find  u  :  Q  — >  ffi,  such  that 

Cu  =  f  in  Q  (156) 

u  =  g  on  P  (157) 

u(0+)  =  it( 0_)  on  H0  (158) 

where,  for  definiteness,  C  is  taken  to  be  time- dependent,  advection- diffusion  operator, 
that  is, 

d 

£  =  —  +  a  •  V  —  kA  (159) 

in  which  n  >  0  is  a  given  constant  and  a  is  a  given  solenoidal  vector  field,  viz. 

V  •  a  =  0  (160) 

Remark 

The  results  to  follow  are  applicable  to  a  wider  class  of  problems  than  the  one  considered 
here.  However,  focusing  on  it  simplifies  the  presentation. 
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Let 

a(w,  u)q  =  -  +  a  •  Vw,  u^j  -  (Vw,  kVu)q 

where  ( u,v)q  =  uv  dQ.  Then 
JO 


a(w,  u)q  =  (w,  Cu)q  =  ( C*w,u)q 


for  all  sufficiently  smooth  w,  u  such  that 


where 


u  =  w  =  0  on  P' 


£*  =  -—-  a  ■  V  -  kA 
dt 


■u  =  u  +  u 
w  =  w  +  w' 


(161) 

(162) 

(163) 

(164) 

(165) 

(166) 


where  we  assume 

v!  =  w'  =  0  on  P'  (167) 

The  variational  equation  corresponding  to  the  initial/boundary- value  problem  is 

{w{T~),u{T~))qt  +  a(w,  u)Q  =  (w,  f)Q  +  (w(0+),  w(0"))no  (168) 

In  (168),  (-,  -)nT  and  (-,-)n0  denote  L2  inner  products  on  O t  and  flo,  respectively.  The 
Euler-Lagrange  equations  corresponding  to  (168)  are  (156)  and  (157).  The  Dirichlet  boundary 
condition,  (157),  is  assumed  satisfied  ab  initio  by  the  trial  solution  u.  Likewise,  w  is  assumed  to 
satisfy  a  homogeneous  Dirichlet  boundary  condition.  Substituting  (165)  and  (166)  into  (168) 
leads  to 


(■ w(T  ),  u(T  ))nT  +  a(w  +  w',u  +  u')q  =  (w  +  w',  /)q  +  (w(0+),m(0  ))n0  (169) 

Note  that  the  integrals  over  Qt  and  Qq  are  unaffected  by  u'  and  w'  because  of  assumption 
(167).  Assuming  w  and  w1  are  linearly  independent,  gives  rise  to  two  subproblems: 

(w(T~),u(T~))nT  +a(w,u)Q  +  a(w,u')Q  =  (w,  f)Q  +  (w(0+),u(0~))no  (170) 

(w(T~),u(T~))nT  +a(w,u)Q  +  (C*w,  u')Q  =  (w,  f)Q  +  (w(0+),u(0~))no  (171) 

a(w',u)Q  +  a(w',u')Q  =  (w',f)Q  (172) 

(w',Cu)q  +  (w',Cu')q  =  (■ w',f)Q  (173) 

In  deriving  (171)  and  (173)  we  used  (162). 

Equation  (173)  is  equivalent  to  the  following  system  of  element-wise  problems  for  u '  \ 


Cu'  =  -(£u-  /) 
u'  =  0 


in  Qe 
on  Pe 


(174) 
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Unfortunately,  these  problems  are  ill-posed  in  the  sense  that  they  are  overspecified.  The 
space-time  Dirchlet  condition  v!  =  0  on  Pe  cannot  be  satisfied  due  to  the  causal  evolutionary 
structure  of  £.  To  circumvent  this  difficulty,  consider  an  elliptic  regularization  of  (174): 


Ceu'£  =  —(£u  —  /) 
u’e  =  0 


in  Qe 
on  Pe 


1,2,...,  ne\ 


(175) 


where 

c‘  =  Tt  +  a-v-KA-si  (176) 

These  problems  are  well-posed  for  all  £  >  0. 

The  solution  of  (175)  can  be  obtained  with  the  aid  of  the  corresponding  Green’s  function 


problems 

£*g'  =  <5  inQe)e_i2 

g'e  =  0  on  Pe  J6  ’ 

(177) 

Thus 

u'£{r0)  =  -[  g'E(r,r0)(£u  -  f)(r)dQr 

(178) 

JQ' 

u'£  =  M'e(£u  -  f ) 

(179) 

where 

r  =  {cc,  t }  G  Q,  r0  =  {*o,  M  €  Q 

(180) 

Remarks 

1.  Cu  —  f  is  the  residual  of  the  resolved  scales. 

2.  The  subgrid  scales  are  driven  by  the  residual  of  the  resolved  scales. 

The  desired  variational  equation  is  obtained  by  taking  the  limit  £  — >  0.  Let 

g'  =  Hmg'  (181) 

£ — ^0 

M'  =  lim  M'  (182) 

£ — ►  O 

Note  that  g'  is  the  causal  Green’s  function  (see,  e.g.,  Stakgold  (1979))  for  the  adjoint 
operator  £* . 

Substituting  u'£  into  (171)  for  u'  and  taking  the  limit  £  — >  0  yields 


(w(T  ),  u{T  ))nT  +  a(w,  u)Q  +  ( C*w ,  M\Cu  -  /))q/  =  (w,  f)Q  +  (w( 0+),  u( 0  ))n0 


(183) 

where 

(£*w,  M'(£u  —  f))Q'  =  -  f  f  £*w(r0)g'(r,r0)(£u- f)(r)dQrdQro  (184) 

JQ'  JQ' 


and 


(185) 


Encyclopedia  of  Computational  Mechanics.  Edited  by  Erwin  Stein,  Rene  de  Borst  and  Thomas  J.R.  Hughes. 
©  2004  John  Wiley  &  Sons,  Ltd. 


49 


Remarks 

1.  Note  that  this  integration  is  over  element  interiors  Q' 

2.  The  effect  of  the  unresolved  scales  on  the  resolved  scales  is  exactly  accounted  for,  up  to 
the  assumption  v!  —  0  on  P1 . 

3.  The  nonlocal  effect  of  the  unresolved  scales  on  the  resolved  scales  is  confined  within 
individual  elements. 

The  variational  equation  (183)  can  be  written  succinctly  as 

(186) 

where: 

B(w,u\  g)  =  (w(T~),u(T~))nT  +  a(w,u)Q  +  (jC*w,M(jCu))qi  (187) 

L{w\  g)  =  (w,f)Q  +  (w{0+),  u(CT))n0  +  (£*«>,  Mf)Q,  (188) 

The  numerical  version  of  (186)  can  be  developed  by  selecting  finite  element  approximations  of 
the  trial  solution  and  weighting  functions,  uh  ss  u  and  wh  ss  w,  respectively,  and  approximating 
the  element  Green’s  function,  g'  ss  g',  viz. 

(189) 

This  amounts  to  applying  Galerkin’s  method  to  the  subgrid-scale  model  obtained  from 
the  variational  multiscale  procedure. 

4-5.  Bubbles  in  space-time 

Now  consider  typical  finite  element  shape  function  expansions  plus  bubbles  (see,  e.g.  Baiocchi, 
Brezzi  and  Franca  (1993),  Brezzi  et  al.  (1992),  Franca  and  Farhat  (1994a, 1994b, 1995),  Franca 
and  Frey  (1992)).  The  set-up  is  as  follows: 


uh 

^nodes 

=  ^  Naua  (likewise  wh) 

(190) 

A=1 

Na 

=  standard  finite  element  shape  functions 

(191) 

UA 

=  nodal  values 

(192) 

^-bubbles 

/ 

u 

=  ^  N'^u  a  (likewise  w') 

(193) 

A= 1 

N'a 

=  bubble  functions 

(194) 

U  a 

=  generalized  coordinates 

(195) 

The  situation  for  linear  finite  element  shape  functions  plus  typical  bubbles  is  illustrated  in 
Figure  21. 
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Substituting  these  functions  into  (171)  and  (173),  and  eliminating  u' a  by  static 
condensation  yields: 

B(wh,uh;g')  =  L(whrg')  (196) 

where 

B(wh,uh;g')  =  (■ wh(T-),uh(T-))nT+a(wh,uh)Q  +  (C*wh,M\Cuh))Q/  (197) 

L(wh;g')  =  (whJ)Q  +  (w(0+),u(0-))no  +  (C*wh,M'f)Q,  (198) 

and 

{C*wh,M'{Cuh-f))Q,  =  -(  [  (£*wh)(ro)g'(r,ro)(Cuh)(r)dQrdQro  (199) 

Jq’  Jq' 

^bubbles 

g'(r,r0)=  ^  KirMN^N^N'M  (200) 

A,B= 1 

This  is  an  approximate  subgrid  scale  model.  Bubbles,  according  to  (200),  generate  an 
approximate  element  Green’s  function,  g'  «  g'.  In  practice,  the  quality  of  the  approximation 
may  be  poor ,  because  standard  polynomial  bubbles  do  not  adequately  represent  the  fine-scale 
structures  which  characterize  u' . 

^.6.  Stabilized  methods 

Stabilized  methods  in  the  present  case  are  generalized  Galerkin  methods  of  the  form 

(wh(T~),  uh (T~))nT  +a(wh ,  uh)Q+(Lwh ,  t(Cuh — /))q'  =  (wh,  f)Q+(w(0+),  u(0~))no  (201) 

where,  typically,  L  is  a  differential  operator  and  r  is  an  algebraic  operator.  Examples  of 
L  are 


L  =  +C 

Galerkin/least-squares  (GLS) 

L  =  +£adv  =  ■§[  +  a,  ■  V 

SUPG 

(202) 

L  =-C* 

Multiscale 

By  comparing  (200)  with  (189),  a  stabilized  method  of  multiscale  type  is  seen  to  be  an 
approximate  subgrid  scale  model  in  which  the  algebraic  operator  r  approximates  the  exact 

integral  operator  M,  i.e. , 

r  =  -M'  «  —M'  (203) 

Equivalently, 

r  •  S(r0  -  r)  =  g'(r,  r0)  «  g'(r,  r0)  (204) 

This  assertion  is  established  by  the  following  calculation: 
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/  /  (-C*wh){r0)g\r,r0)(Cuh  -  f)(r)dQrdQro 

Jq '  JQ’ 

[  (  (- C*wh)(r0)T6(i'0  -  r)(£uh  -  f){r)dQrdQr 

JQ'  JQ' 

[  (~£*wh)(r)T(£uh  -  f)(r)dQr 

JQ' 


(205) 


£6.1.  Formulas  for  r  Formulas  for  r  may  be  derived  from  (204).  For  example,  assume  it  is 
sufficient  to  approximate  r  as  a  constant  over  each  element.  Then 


/  /  rS(r0  -  r)dQrdQro  =  g'(r,  r0)dQrdQro 

JQe  JQ e  JQ"  JQ" 

g '{r,  r0)dQrdQro 


I  Qe  JQe 

I  I 

)Q‘  JQ" 


meas (Qe)  JQ e  7Q, 


g'(r,  r0)dQrdQr 


(206) 

(207) 


4-6.2.  Example:  First-order  ordinary  differential  equation  in  time  To  demonstrate  the  use  of 
formula  (207),  consider  the  simple  example  of  a  first-order  ordinary  differential  equation  in 
time.  (In  Section  3  the  exact  value  of  r  for  the  steady  advection-diffusion  operator  in  space 
was  determined.)  Let 

£  =  |  (»8) 

£e  =  £-e^  (209) 

The  initial-value  problem  is:  Given  /  :  Q  — >  ffi,  u(0“)  G  M,  find  u  :  Q  — >  ffi  such  that 


£u  =  f  in  Q  =  [0,  T] 
u(0+)  =  u(0") 

The  element  Green’s  function  problems  for  the  regularized  equation  are 

«  =  <5 

ge  =  0 

where 

d_ 

dt,  “  dt2 )  dt  dt2 


in  Qe 
on  dQe 


e  —  1,2,...,  ne  i 


d2\ 


d  d2 


(210) 

(211) 


(212) 


(213) 


The  element  Green’s  function  is  sketched  in  Figure  28  as  a  function  of  a  =  At/(2e),  where 
At  =  meas(<5e)  =  element  “length”.  Note  that  the  limit  a  — »  oo  (equivalently  £  — >  0),  gives 
rise  to  the  causal  Green’s  function  for  £*. 
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oo 


a  >  1 


a  <  1 


0 


Figure  28.  Green’s  function  for  The  limit  a  — >  oo  is  the  causal  Green’s  function  for  C* . 


Let 


Then 


Ta 


1 

meas(Qe) 


g  (t,t0)dQtdQt0 


At 

~2 


^coth  a  — 


t  =  lim  Ta 
06— >00 


At 

~2 


(214) 


(215) 


Remark 

This  expression  for  r  is  known  to  be  the  correct  one  for  the  case  at  hand  from  the 
mathematical  theory  of  stabilized  methods.  It  is  clear  that  this  value  of  r  is  as  good  as  using 
the  exact  Green’s  function  representation  as  long  as  the  trial  functions  are  piecewise  constant 
with  respect  to  time  on  element  subdomains. 


4-7.  Summary 

The  development  of  a  class  of  subgrid  scale  models  for  time-dependent  problems  was  presented. 
The  corresponding  steady  case  was  described  previously  in  Section  3.  It  was  shown  that  bubble 
functions  give  rise  to  approximate  element  Green’s  functions  in  the  subgrid-scale  models,  albeit 
typically  not  good  approximations.  Likewise,  stabilized  methods  were  identified  as  approximate 
subgrid-scale  models.  Formulas  for  r  emerged  from  this  identification  and  this  opens  the  way 
to  improved  representations  of  r  in  stabilized  methods  for  equations  of  evolution.  Stabilized 


Encyclopedia  of  Computational  Mechanics.  Edited  by  Erwin  Stein,  Rene  de  Borst  and  Thomas  J.R.  Hughes. 
©  2004  John  Wiley  &  Sons,  Ltd. 


53 


methods  have  been  shown  to  be  “good”  numerical  methods  for  problems  in  which  the  classical 
Galerkin  finite  element  method  fails  (see,  e.g.  Barbosa  and  Hughes  (1992),  Codina  (1998,2000), 
Franca  and  Do  Carmo  (1989),  Franca  and  Farhat  (1995),  Franca  and  Frey  (1992),  Franca, 
Frey  and  Hughes  (1992),  Franca  and  Hughes  (1988),  Franca  and  Hughes  (1993),  Franca  et  al. 
(1988),  Franca,  Hughes  and  Stenberg  (1993),  Franca  and  Valentin  (2000),  Hughes  and  Brezzi 
(1989),  Hughes  and  Franca  (1987),  Hughes,  Franca  and  Balestra  (1986),  Hughes,  Franca  and 
Hulbert  (1989),  Hughes,  Hauke  and  Jansen  (1994)  in  which  error  estimates,  stability  results, 
and  verification  problems  are  presented) .  The  results  described  in  this  and  the  previous  section 
provide  an  alternative  perspective  of  stabilized  methods. 

Starting  with  partial  differential  equation  plus  boundary  and  initial  conditions,  application 
of  the  variational  multiscale  method,  then  the  Galerkin’s  method,  results  in  a  numerical  method 
that  has  the  form  of  a  stabilized  method.  In  addition,  it  provides  detailed  expressions  that 
enable  one  to  derive  more  accurate  stabilized  methods. 


5.  Stabilized  Methods  for  Aclvective-Diffusive  Equations 

Stabilized  methods  for  steady  and  unsteady  advective-diffusive  systems  are  described.  Fairly 
general  boundary  conditions,  leading  to  well-posed  variational  problems,  are  considered.  The 
boundary-value  problems  are  specialized  to  the  hyperbolic  case  for  completeness.  Galerkin, 
SUPG,  Galerkin/least-squares,  and  multiscale  stabilized  finite  element  methods  are  contrasted. 
An  a  priori  error  analysis  of  Galerkin  least-squares  is  presented.  The  developments  for  the 
steady  case  are  generalized  to  the  unsteady  case  by  way  of  space-time  formulations  employing 
the  discontinuous  Galerkin  method  with  respect  to  time.  Symmetric  advective-diffusive  systems 
are  also  considered.  The  set-up  so  closely  follows  the  scalar  case  that  completely  analogous 
results  are  obtained.  In  order  to  fully  comprehend  these  developments,  the  reader  is  urged  to 
first  carefully  study  the  scalar  case,  as  the  system  case  is  presented  in  virtually  equation-for- 
equation  form  with  little  amplification. 

5.1.  Scalar  steady  advection- diffusion  equation 

5.1.1.  Preliminaries  Let  O  be  an  open,  bounded  region  in  ffid,  where  d  is  the  number  of  space 
dimensions.  The  boundary  of  is  denoted  by  T  and  is  assumed  smooth.  The  unit  outward 
normal  vector  to  T  is  denoted  by  n  =  (ni,n2,  ■  ■  ■  ,n<i )•  Let  a  denote  the  given  flow  velocity, 
assumed  solenoidal,  that  is,  V  •  a  =  0.  The  following  notations  prove  useful: 


Let  {r-,r+}  and  {Tg,Th} 

an  —  n  a 

an  (an  J-  \an 

a~  =  (an  -  |  an 

be  partitions  of  T,  where 

l)/2 

l)/2 

(216) 

(217) 

(218) 

T-  = 

{x  G  T  |  an(x)  <  0} 

(inflow  boundary) 

(219) 

r+  = 

r-r- 

(outflow  boundary) 

(220) 

The  following  subsets  are  also  needed  (see  Figure  29): 
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Figure  29.  Illustration  of  boundary  partitions. 

r9  =  (221) 

it  =  r^nr±  (222) 

Let  k  =  const.  >  0  denote  the  diffusivity.  Various  fluxes  are  employed  in  the  sequel: 


0) 

=  —au 

(advective  flux) 

(223) 

(w) 

=  k\7u 

(diffusive  flux) 

(224) 

(T 

=  cra  +crd 

(total  flux) 

(225) 

=  n  ■  (Ta 

(226) 

_  d 
an 

=  n  ■  crd 

(227) 

=  n  ■  (T 

(228) 

Let  D  denote  a  domain  (e.g.,  Cl,  T,  etc.).  The  L,2(D)  inner  product  and  norm  are  denoted 
by  (•,  -)d  and  ||  ■  ||_d,  respectively. 

5.1.2.  Problem  statement  The  problem  consists  of  finding  u  =  u{ x)  \/x  £  Cl,  such  that 


s 

III 

<1 

II 

in  Cl 

(229) 

u  =  g 

on  Tg 

(230) 

a~u  +  ad(u)  =  h 

on  Tft 

(231) 
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where  /  :  Q  — >  ffi,  g  :  Tg  — ->  ffi,  and  h  :  — >  ffi  are  prescribed  data.  (229)  is  an  elliptic 

equation.  The  boundary  condition  (231)  can  be  better  understood  by  letting 


h  = 


h  on  Th 
h+  on  T+ 


Thus,  (231)  may  be  written  in  the  equivalent  form: 

cr„(«)  =  h~  on  T^  (total  flux  b.c.) 

<Tn(u)  =  on  (diffusive  flux  b.c.) 


(232) 


(233) 

(234) 


5.1.3.  Variational  formulation  The  variational  form  of  the  boundary- value  problem  is  stated 


in  terms  of  the  following  function  spaces: 

5  =  {u€H1(Q)\u  =  gonTg}  (235) 

V  =  {w€H1(Q)\w  =  0onTg}  (236) 

The  objective  is  to  find  u  €  S  such  that  Vic  G  V 

B(w,u)  =  L(w)  (237) 

where 

B(w,u)  =  (Vw,(r(u))n  +  (w,a+u)rh  (238) 

L(w)  =  (w,f)n  +  (w,h)rh  (239) 


The  formal  consistency  of  (237)  with  the  strong  form  of  the  boundary- value  problem,  that 
is,  (229)-(231),  may  be  verified  as  follows: 


0  =  B(w,u)  —  L(w) 

=  ~{w,  V  •  <r{u))n  +  (w,  on(u))rh  +  (w,  a+u)rh  -  (w,  f) n  -  (ic,  h)Th 
=  -(ic,V  -tr(u)  +  f)a  +  (w,-a~u  +  a^(u)  -  h)Th  (240) 

Stability ,  or  coercivity ,  is  established  as  follows 


B(w,  w) 


(Vic,  —aw  +  kVw)q  +  (w,  a+w)rh 
-^(w,anw)rh  +  «||Vw||n  +  (w,a+w)Th 


cll  Vic||j 


2 

r„ 


Vic  €  V 


(241) 


For  future  reference,  define: 

|||ic|||2  =  B(w,w)  (242) 

Finally,  the  global  conservation  of  flux  is  investigated.  Consider  the  case  in  which  Tg  =  0. 
Set  w  =  1  in  (237): 


0 


(243) 
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which  may  be  written  equivalently  as 

0  =  /  h~dT+  [  fcin+  f  (~anu+  h+)  dT,  (244) 

Jr-  Jn  J r+ 

This  confirms  the  conservation  property  for  the  case  assumed.  If  Tg  ^  0,  “consistent”  fluxes 
on  Tg  may  be  defined  via  a  mixed  variational  formulation  which  automatically  attains  global 
conservation.  See  Hughes  (1987)  p.  107,  Hughes  et  al.  (2000)  and  Hughes  et  al.  (1987)  for 
background. 


5.1.4-  Hyperbolic  case  In  the  absence  of  diffusion,  a  boundary  condition  on  the  outflow 
boundary  cannot  be  specified.  The  equations  of  the  boundary- value  problem  are 


Cu  =  -V  •  (Ta{u)  =  f  in  O  (245) 

u  =  g  on  T~  (246) 

<(u)  =  h~  on  T^  (247) 

The  variational  operators  are  defined  as 

B(w,u )  =  (Vu;,  cra(u))n  +  (w,a+u)r  (248) 

L{w)  =  (w,f)n  +  (w,h~)r-  (249) 

Consistency,  stability  and  conservation  are  established  as  follows: 


Consistency 


0  =  B(w,  u)  —  L(w) 

=  -(w,  V  •  cra(u))n  +  (w,  -anu)r  +  (w,  a+u)r  -  (w,  f)n  -  (w,  h~)r- 

=  -(w,  V  •  cra(u)  +  f) n  +  (w,  —a~u  -  h~)r- 

Stability 


Conservation  (T 


Equivalently, 


B(w,w)  =  (Vtt;,  —aw)n  +  (w,a+w)r 
=  -i(u>,  anw)r  +  (w,a+w)r 


0) 


an  2  w 


\/w  £  V 


0  =  B(l,u)-L(l) 


/  a„  dT  f  dfl  —  h  dT 
/  r  Jn  J rr 


0  =  h  dT+  fdil+  —anu  dT 
Jr  Jn  J r+ 


(250) 


(251) 


(252) 


(253) 
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5.1.5.  Finite  element  formulations  Consider  a  partition  of  H  into  finite  elements.  Let  fle  be 
the  interior  of  the  eth  element,  let  Te  be  its  boundary,  and 


e 

(element  interiors) 

(254) 

Ure\r 

(element  interfaces) 

(255) 

e 


Let  Sh  C  S ,  Vh  C  V  be  finite  element  spaces  consisting  of  continuous  piecewise  polynomials 
of  order  k.  As  a  point  of  departure,  consider  the  classical  Galerkin  method: 

Find  uh  €  Sh  such  that  \/wh  €  Vh 

B(wh,uh)  =  L(wh)  (256) 

Remark 

The  element  Peclet  number  is  defined  by  a  =  h\a\/[2n).  The  entire  range  of  a,  that  is, 
0  <  a  <  oo,  is  important.  The  advection  dominated  case  (i.e.,  a  large)  is  viewed  as  “hard”. 
The  Galerkin  method  possesses  poor  stability  properties  for  this  case.  Spurious  oscillations  are 
generated  by  unresolved  internal  and  boundary  layers.  The  reason  B(w,u)  is  not  capable  of 
suppressing  the  spurious  oscillations  can  also  be  easily  inferred  from  its  stability  bound  (251). 
This  bound  indicates  that  the  bilinear  form  of  the  Galerkin  method  does  not  exercise  any 
control  over  the  first  derivatives  of  the  solution.  As  a  result,  the  norm  of  the  gradient  of  the 
solution  can  grow. 

Methods  with  improved  stability  properties  are  given  below: 

SUPG 

BSUPG{wh,Uh)  =  LSUPG(wh ) 

BSuPG{wh,  uh)  =  B(wh,uh)  +  ( ra  ■  Vwh,  Cuh)n' 

Lsupc(wh)  =  L(wh)  +  (ra  •  Vwh,  /)q/ 

Galerkin/least-squares 

BGLS(wh,uh)  =  LGLS(wh) 

BGLs{wh,  Uh)  =  B(wh,  uh)  +  ( rCwh ,  Cuh)n’ 

LGLs{wh )  =  L(wh)  +  (rCwh,  f)n> 

Multiscale 

BMS(wh,uh)  =  LMs{wh ) 

BMS(wh,uh)  =  B(wh,uh)-(T£*wh,jCuh)n> 

LMS(wh)  =  L(wh)-(TC*whJ)n > 

C*(wh)  =  -a-Vwh—K  Awh  (A  is  the  Laplace  operator) 

Remarks 

1.  In  the  hyperbolic  case,  or  for  piecewise  linear  elements  in  the  general  case,  SUPG, 
Galerkin/least-squares,  and  multiscale  become  identical. 

2.  SUPG,  Galerkin/least-squares,  and  multiscale  stabilized  methods  are  residual 
methods.  That  is,  (257),  (260),  and  (263)  are  satisfied  if  uh  is  replaced  by  u,  the 
exact  solution  of  the  boundary- value  problem. 


(257) 

(258) 

(259) 

(260) 
(261) 
(262) 

(263) 

(264) 

(265) 

(266) 
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Figure  30.  Definition  of  C(a). 


5.1.6.  Error  analysis  The  SUPG  method  was  originally  analyzed  in  Johnson,  Navert 
and  Pitkaranta  (1984)  and  Navert  (1982).  In  this  section  an  a  priori  error  analysis  of 
Galer kin/least-squares  is  performed.  The  multiscale  method  may  be  analyzed  using  similar 
techniques. 

Let  e  =  uh  —  u  denote  the  error  in  the  finite  element  solution.  By  Remark  2,  above, 

BGLS{wh,e)  =  0,  VwheVh  (267) 

This  is  referred  to  as  the  consistency  condition  for  Galer  kin/least-squares. 

Let 

(268) 

By  (261)  and  (268), 

BGLs{wh,wh)=\\\wh\\\2GLS,  Wwh  e  Vh  (269) 

This  is  the  stability  condition  for  Galerkin/least-squares. 


2 

GLS 


-2Cwh 


Remarks 

1.  Stability  is  less  straightforward  for  SUPG  and  multiscale.  One  needs  to  invoke  an 
“inverse  estimate”  and  specific  properties  of  r.  These  assumptions  are  seen  to  be 
unnecessary  for  establishing  the  stability  of  Galerkin/least-squares.  However,  they 
resurface  in  the  convergence  analysis. 
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2.  A  term  of  the  form  ||uU||q  can  be  added  to  (268)  by  employing  a  change  of  variables. 
See  Johnson,  Navert  and  Pitkaranta  (1984)  and  Navert  (1982)  for  further  discussion. 

3.  For  the  related  residual-free  bubble  approach  (that  coincides  with  SUPG,  with  a  specific 
choice  of  r  for  linear  elements,  but  is  slightly  different  for  quadratic  and  higher-order 
elements)  error  estimates  were  proved  in  Brezzi  et  al.  (1999)  and  Brezzi,  Marini  and  Siili 
(2000) .  Additional  results,  including  local  error  bounds  were  proved  in  Sangalli  (2000) . 
See  also  Asensio,  Russo,  and  Sangalli  (2004)  for  additional  comments  on  residual-free 
bubbles  for  quadratic  elements. 

4.  The  quasi-optimality  of  SUPG  methods  with  respect  to  the  suitable  problem-dependent 
“ideal”  norms  was  analyzed  in  depth  by  Sangalli  (2004). 

5.  SUPG  and  the  related  residual-free  bubble  methods  were  also  analyzed  from  the  point 
of  view  of  a  posteriori  error  estimates  in  Russo  (1996b),  Verfiirth  (1998),  Papastavrou 
and  Verfiirth  (2000),  Sangalli  (2001),  and  references  therein. 


Let  uh  £  Vh  denote  an  interpolant  of  u.  The  interpolation  error  is  denoted  by  rj  =  u"  —  u. 
Thus,  e  =  eh  +  rj,  where  eh  £  Vh. 

It  is  assumed  that  r  possesses  the  following  properties: 

a  large  (270) 

a  small  (271) 

A  specific  choice  of  r  satisfying  these  properties  is  given  by 


where  (((a)  is  illustrated  in  Figure  30.  (See  Hughes,  Mallet  and  Mizukami  (1986),  Appendix  I, 
for  some  other  possibilities.) 

For  sufficiently  smooth  u,  standard  interpolation  theory  (see,  e.g.,  Ciarlet  (1978))  and  the 
above  asymptotic  properties  of  r  lead  to  the  following  interpolation  estimate : 


rh 


r2Cr] 


<  cuh 


(273) 


J  2k +  1,  a  large 
(  2k,  a  small 


(274) 


where  cu  is  a  function  of  u.  The  notation  cu  is  used  subsequently,  it  being  understood  that  in 
each  instance  its  value  may  change  by  a  multiplicative  constant. 

An  inverse  estimate  also  needs  to  be  introduced.  The  appropriate  form  in  the  present 
circumstances  is 

||Aioh||n,  <c/»_1||Vioh||n  \/wh£Vh  (275) 

where  c  is  a  nondimensional  constant.  (See  Ciarlet  (1978),  pp.  140-146,  for  results  of  this 
kind.) 
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Theorem  5.1.  Assume  the  consistency  condition  (267),  stability  condition  (269),  and 
interpolation  estimate  (273)  hold.  Assume  the  slope  m  in  the  definition  of  £(a)  satisfies 
m  <  4 /c,  where  c  is  the  constant  in  the  inverse  estimate  (275).  Then  the  error  for  the 
Galerkin/least-squares  method  is 

IINIIgls  —  °uh21  (276) 

PROOF. 

First,  estimate  eh: 


I GLS 


=  BGLS(e\eh) 

=  BGLS(eh,e-  ??) 
=  - BGLS(eh,r] ) 

<  \BGLS{eh,ri)\ 


(stability) 

(consistency) 


| -(a  •  Veh,rj)n  +  n{\7eh,  V??)n 

+(eh,  a+??)rfe  +  {rCeh,  Crj)w  |  (definition  of  BGLS(-,  •)) 
\-(Ceh,rj)n'  -  n(Aeh,  77)^ 

+«(Veh,  V??)n  +  (eh,  a+rj) Th 
+  {rCeh,  £??)□/ 1 


< 


r2Ce 
2 


n 


T2  Aeh 
2 


1 

'I 
1 

'I 


K  II V77H 


r  277 
2 


n 


|a„.  1 2 


-2Ceh 


2 

1  1 1 

+ 

\an\2V 

rh 

2 

n 

2 


t  2  Trj 


To  proceed  further,  invoke  the  bound  on  ?n; 


nh 

2R 

4a 

h2 

<  — 


C(«) 


C(a) 


(277) 


(278) 


Combining  (278)  with  the  inverse  estimate  yields 


7-2  Aeh 


Employing  the  result  in  (277)  leads  to 


<2 
I  GLS  ~ 


<  k  Ve'1 


:l|Vr/||( 


|an|2?? 


rh 


t  2  Trj 


(279) 


(280) 
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Therefore,  by  the  interpolation  estimate, 

llle^lll2  <ch21 

111°  III QLjS  —  L'u  0 

(281) 

Likewise, 

\MfGLS<Cuh21, 

(282) 

and  so,  by  the  triangle  inequality, 

lllclll2  <  c  h21 

lllelllGLS  —  <uU 

(283) 

This  completes  the  proof  of  the  theorem.  □ 

Remark 

The  same  rates  of  convergence  can  be  proved  for  SUPG  and  multiscale  stabilized  methods. 


5.2.  Scalar  unsteady  advection- diffusion  equation:  Space-time  formulation 

The  initial/boundary-value  problem  consists  of  finding  u(x,t),  Vx  £  O,  Vt  £  [0,T],  such  that 


Ctu  =  ii  +  Cu=f  in  S2x]0,  T[  (284) 

w(x,0)  =  uo(£c)  \/x  £  O  (285) 

u  =  g  onTsx]0,T[  (286) 

— afu  +  a^(u)  =  h  onThx]0,T[  (287) 


where  ii  =  du/dt,  and  uq  :  SI  — *  ffi,  /  :  flx]0,T[— *  R,  g  :  Tgx]0,T[— >  R,  and  h  :  Tft,x]0,T[— >  R 
are  prescribed  data. 

The  procedures  to  be  described  are  based  upon  the  discontinuous  Galerkin  method 
in  time.  (See  Johnson  (1987),  and  references  therein,  for  a  description  of  the  discontinuous 
Galerkin  method.)  Space-time  (i.e.  f2x]0,T[)  is  divided  into  time  slabs  fix]f„, [,  where 
0  =  to  <  ti  <  ■  ■  ■  <  tN  =  T.  Each  time  slab  is  discretized  by  space-time  finite  elements.  The 
finite  element  spaces  consist  of  piecewise  polynomials  of  order  k  in  x  and  t,  continuous  within 
each  slab,  but  discontinuous  across  time  slabs.  Again,  as  a  point  of  departure,  the  Galerkin 
method  will  be  presented  first: 


B(wh,uh)n 
B(wh,  Uh)n 


L(wh)n 


L(wh)n  ,  n  =  0, 1, . . . ,  IV  —  1 

rtn  +  1 

/  (~wh,  uh)n  +  B(wh,uh))  At  +  (wh  (t~+1),  uh{t~+1))n 


L(w  )dt+(w  (t+),u  (tn))a 


(288) 

(289) 

(290) 

(291) 


where 


Uh(tn)  =  uh(x^n)  (292) 

uh(fo)  =  u0(x)  (293) 


Remark 

Continuity  of  the  solution  across  time  slabs  is  seen  to  be  weakly  enforced. 
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Generalization  of  SUPG,  Galer kin/least-squares  and  multiscale  proceeds  analogously  to  the 
steady  case: 


SUPG 


BsUPG(wh  ,Uh)n  = 

LsUPG{wh)n  , 

3 

II 

o 

H 

1 

1 — 1 

(294) 

BsUPG(wh  ,Uh)n  = 

B{ 

Pl'Tl- 

Wh  ,  Uh)n  +  J 

(r  (wh  +  a  •  Vioil)  ,  Ctuh)Q,  dt 

(295) 

LsUPG{wh)n  = 

rtn+1 

L(wh)n  +  J  (l 

r(wh  +  a-Vwh),f)a,  dt 

(296) 

Galerkin/least-squares 

BGLs(wh,Uh)n 

= 

LGLs{wh)n  , 

n  =  0,1, . . .  ,N  —  1 

(297) 

BGLs{wh,Uh)n 

= 

B(wh,Uh)n  + 

rtn  +  1 

J  (TCtwh,Ctuh)n,  dt 

(298) 

LGLs(wh)n 

= 

L(wh)n  +  J 

r+1 

(TCtwh,f)n,  dt 

(299) 

Multiscale 

BMs(wh,  Uh)n 

= 

LMs{wh)n  > 

n  =  0,1, . . .  ,N  —  1 

(300) 

BMS(wh,  Uh)n 

= 

B{wh,uh)n- 

rtn+1 

J  (TC*twh,Ctuh)n,  dt 

(301) 

LMs{wh)n 

= 

L(wh)n  -  j  “ 

(■ rC;wh,f)a,dt 

(302) 

r* 

= 

-wh  +  C*wh 

(303) 

Remarks 


1.  In  the  unsteady  case,  h  represents  a  space-time  mesh  parameter. 

2.  The  issue  of  the  time  integration  method  is  obviated  by  the  choice  of  space-time 
interpolation.  Unconditional  stability  (i.e.,  stability  for  all  time  steps,  At  =  tn+ 1  —  fn  > 
0)  is  achieved  in  all  cases.  On  each  time  slab  a  system  of  linear  algebraic  equations 
needs  to  be  solved. 


Let 


JV-l 


I  I2  =  5  £  l|[«’‘(Wllln  +  \  (l|«’'*<r-)lln  +  IK(0+)lln) 


n—  1 
cT 


J  d  t 


where  [wh(tn)]  =  wh(t^)  —  wh(tn).  It  is  a  simple  exercise  to  show  that 


N- 1 


N—  1 


Y  B(wh,wh)n  -  ^(w/l(t+),w/l(t„))n  =  i  Wh 


n— 0 


(304) 


(305) 
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where  the  constraint  uih(t0  )  =  wh{ 0  )  =  0  has  been  used. 

Likewise, 

JV-l  N- 1 

E  BGLs(wh,Wh)n  -  E  =  |  Wk  |  2GLS  (306) 

71=0  71=  1 

where 


w 


2 

GLS 


=  u? 


N-l 

t+e 

n— 0  ' 


^n+Jl  n  .  ,  |i  2 


T2CtW 


dt 


(307) 


The  following  error  estimate,  analogous  to  the  steady  case,  can  be  established  for  the  space- 
time  Galerkin/least-squares  method: 


e 


GLS 


<  c„h21 


(308) 


Remark 

The  hypotheses  necessary  to  prove  (308)  are  virtually  identical  to  those  for  the  steady  case. 


5.3.  Symmetric  advective- diffusive  systems 

The  previous  developments  for  the  scalar  advection-diffusion  equation  may  be  generalized  to 
symmetric  advective-diffusive  systems.  The  equations  are  (see  also  Hughes,  Franca  and  Mallet 
(1987),  Hughes  and  Mallet  (1986)): 


CtV 

CV 

V 


k 


A  ■  VV 


A  0V,t+£V  =  f 
A  ■  VV  -  V  •  KSJV 


• ,  vn)T 

A-i,  A.2, 

•  •  i  Ad 

'  kn 

...  k 

.  kd i 

...  k 

~t  ~  ~  av  ~  av 

A  VV  =  At VV,j  =  A1—  +  ...  +  Ad— 

OX 1  OXd 


(309) 

(310) 

(311) 

(312) 

(313) 


(314) 


in  which  Aq  is  an  to  x  m  symmetric,  positive-definite  matrix;  A,  is  an  in  x  to  symmetric, 
matrix,  1  <  *  <  d,  A  is  a  solenoidal  matrix,  that  is,  A.,^  =  0;  and  K  is  an  ( m  ■  d)  x  {in  ■  d) 
constant,  symmetric,  positive-definite  matrix.  (The  case  in  which  K  is  positive- semidefinite  is 
more  important  in  practice,  but  complicates  the  specification  of  boundary  conditions.) 
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Corresponding  to  the  developments  for  the  scalar  case,  the  following  results  are  obtained: 


A 

■**-n 

=  niAi 

(315) 

=  2  +  An  ^ 

(316) 

K 

(317) 

U(V) 

=  A0y 

(temporal  flux) 

(318) 

F?(V) 

=  -AiV 

(advective  flux) 

(319) 

Ff(V) 

=  -KijVj 

(diffusive  flux) 

(320) 

F 

=  Ff(V)  +  Ff(V) 

(total  flux) 

(321) 

K 

=  ruF* 

(322) 

K 

=  ruFi 

(323) 

F 

1  n 

=  riiFi 

(324) 

For  simplicity,  assume  that  for  x  G  T,  An(x)  is  either  positive-  or  negative-definite.  This  will 
allow  a  concise  statement  of  boundary  conditions  analogous  to  the  scalar  case.  For  situations 
in  which  An  is  indefinite,  boundary  condition  specification  is  more  complex,  necessitating 


component  by  component  specification.  Let 

=  {xer  An(x)  < o  j  (325) 

r+  =  r  -  r~  (326) 

rf  =  rgflr*  (327) 

r«  =  r«nr±  (328) 

5.3.1.  Boundary-value  problem 

CV  =  — V  •  F  =  -Fiti  =  T  on  Q  (329) 

V  =  Q  on  Tg  (330) 

-A~V  +  F*(V)  =  n  onVH  (331) 

(331)  is  equivalent  to 

Fn(V)  =  on  (total  flux  b.c.)  (332) 

F*(V)  =  7~L+  on  (diffusive  flux  b.c.)  (333) 

variational  formulation 

s  =  {VcH\n)m  \V  =  g  on  Tg }  (334) 

V  =  {W  G  =  0  on  rH  }  (335) 

B(W,V)  =  (\7W,F(V))n  +  (W,A+V)rn  (336) 

L(W)  =  (W,F)n  +  (W,H)  rn  (337) 


Encyclopedia  of  Computational  Mechanics.  Edited  by  Erwin  Stein,  Rene  de  Borst  and  Thomas  J.R.  Hughes. 
©  2004  John  Wiley  &  Sons,  Ltd. 


65 


0  =  B(W,  V)  —  L(W) 

-  ~(W,  V  •  F(V)  +  F)a  +  (W, -A~V  +  Fdn{V)  -  H) rn 

(formal  consistency) 


(338) 


B(W,W)  = 


K 


VW 


w 


MW  £  V 


lr„ 

(stability) 


\W\\\2  =  B(W,W) 


(339) 

(340) 


0  =  B(1,V)-L(1) 

=  -(y  H~dT  +  Jfdn  +  J  {—AnV  +  H+)  dr^j 

(conservation  for  Tg  =  0)  (341) 

hyperbolic  case 

— V  •  Fa(V)  =  T  on  Cl  (342) 

V  =  g-  on  Tg-  (343) 

K(V)  =  n-  on  IV  (344) 

B(W,V)  =  (M7W,Fa(V))n  +  (W,A+V)  r+  (345) 

L(W)  =  (W,F)n  +  (W,H-)rn_  (346) 


0  =  B(W,  V)  —  L(W) 

=  -(W,S/-Fa(V)  +  F)n  +  (W,~A~V-n~) rH_ 

(formal  consistency)  (347) 


B(W,W) 


MW  e  v 

(stability) 


(348) 


0  =  /  H~dT  +  [  Fd fi  +  [  -AnV  dr 

Jr-  Jn  J r+ 

(conservation  for  Tg-  =  0)  (349) 
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finite  element  formulations 


BsUPG(Wh,Vh )  = 

Lsupg^W11) 

(350) 

BSUPG(W\Vh)  = 

B(Wh,Vh)+  (rA-VWh,CVh^/ 

(351) 

Lsupg{W1)  = 

L(Wh)  +  (tA- 

(352) 

BGLs(Wh,Vh )  = 

LGLS(Wh ) 

(353) 

BGLS(W'\Vh)  = 

B(wh,vh )  +  (rcwh,cvh') 

(354) 

LGLS(Wh )  = 

l(w,1)  +  (tCw'\b ) 

(355) 

BMS(Wh,Vh)  = 

Lms(W'1) 

(356) 

BMS(W\Vh)  = 

B(Wh,Vh )  -  (r£*W,l,£V'*) 

(357) 

LMS(Wh )  = 

L(Wh)  -  (r£*Wh 

(358) 

C*Wh  = 

-AiW$  -  kijW% 

(359) 

GLS-norm  and  error  estimate 

Wh  2  =  BGLS(Wh,Wh) 

GLS 

=  Wh  +  T?CWh  (360) 

O' 

m\\2GLS  <  Cyh21  (361) 


5.3.2.  Initial/boundary-value  problem 


CtV 

=  C7(V)  +  £y  =  T 

in  fix]0,T[ 

(362) 

U{V{x,  0)) 

=  U(Vo(x)) 

Va:  €  0 

(363) 

V 

=  g 

on  rax]0, T[ 

(364) 

A^v  +  F^y) 

=  H 

on  T-h 

(365) 

finite  element  formulations 

B(W}\Vh)n  =  J  ^  ^-Wh,U(Vh)^  +  B(Wh,Vh)jdt 

+  (wh(t-+1),U(yh(t-+S)n  (366) 

L{Wh)n  =  J*n+1  L(Wh)dt+  (wh(t+),U(Vh(t-))^  (367) 
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N- 1  N—l 

wh  |  =  ^2B(Wh,Wh)n-  B(Wh(t+),U(Vh(t-))n 

n—0  n= 1 

1  iV_1  2 

=  2E|4l^(Ul|n 

n=  1 

+  E  Wh(T~)  2+  a|w(0+)  2)+  f  Wh  2clt  (368) 

z  \  ci  ci )  in 


BsupG(Wh,Vh)n  =  LSUPG(Wh)n  ,  n  =  0, 1, . . . ,  IV  —  1  (369) 

BSUPG(Wh,Vh)n  =  B(Wh,Vh)n  +  £  "+1  [r[A0Wh  +  A-\7Wh)  ,  <$70) 

isiAPG(Wh)n  =  L(Wh)n  + Jt"+1  (r(A0Wh +  A-VWh)  (371) 

BGLs(W'\Vh)n  =  LGLS(Wh)n  ,  n  =  0, 1, . . . ,  iV  —  1  (372) 

BGLs(W}\Vh)n  =  B(Wh,Vh)n  +  j^+1  (rCtWh,CtVh^  dt  (373) 

LGLs{Wh)n  =  L(Wh)n  +  J*n+1  (rCtWh,f)^dt  (374) 

BMs(Wh,Vh)n  =  LMS(Wh)n,  n  =  0, 1, . . . ,  IV  —  1  (375) 

BMs(Wh,Vh)n  =  B(Wh,Vh)n-jt  ^  (rC*tWh,CtVh)ndt  (376) 

LMS(Wh)n  =  L(Wh)n-J^+1  (rC*tWh,J^^dt  (377) 


GLS-norm  and  error  estimate 

Likewise, 

iV-l  iV-l 

I  ^  I  gls  =  E  bgls(W\  Wh)n  -  Y,  B(Wh(t+),  U(Wh(t~)))n 

n= 0  n=l 

2  1  /»£n.  +  l  2 

=  I  wM  +  V  /  T?CtWh  dt  (378) 

I  £  I  gls  <  (379) 

Remarks 

1.  Stabilized  methods  have  been  used  widely  in  engineering  applications.  There  is  a  very 
large  literature  on  mathematical  and  practical  aspects.  Experience  has  indicated  that 
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the  SUPG  and  multiscale  variants  are  superior  to  Galerkin/least-squares  in  practical 
applications  (see,  e.g.,  Bochev  and  Gunzburger  (2003)).  A  recent  evaluation  of  SUPG, 
and  comparison  with  finite  volume  and  discontinuous  Galerkin  methods,  is  presented 
in  Venkatakrishnan  et  al.  (2003). 

2.  There  are  interesting  examples  of  non-residual  based  stabilized  methods.  See,  for 
example,  Codina  and  Blasco  (2000a,  2000b),  and  Bochev  and  Dohrman  (2004). 


6.  Turbulence 

In  Sections  2-5,  the  variational  multiscale  method  was  described,  and  its  connection  with 
stabilized  finite  element  methods  was  established.  The  idea  of  a  variational  framework  for 
subgrid-scale  modeling  was  also  discussed.  In  this  Section,  the  application  of  the  multiscale 
method  to  the  incompressible,  Navier-Stokes  equations  is  considered.  The  objective  is  a 
satisfactory  interpretation/generalization  of  the  Large  Eddy  Simulation  (LES)  concept  within  a 
variational  formulation  of  the  Navier-Stokes  equations.  This  requires  dealing  with  nonlinearities 
as  well  as  issues  of  turbulence  modeling. 

The  classical  LES  formulation  of  the  incompressible  Navier-Stokes  equations  is  reviewed 
first.  As  points  of  reference,  filtering,  the  subgrid-scale  stress  and  the  Smagorinsky  model 
are  discussed.  The  estimation  of  Smagorinsky  parameters  by  way  of  the  approach  due  to 
Lilly  (1966,1967,1988)  is  also  recalled.  The  shortcomings  of  the  classical  approach,  noted 
previously  in  the  literature  are  summarized.  A  fundamental  problem  of  the  classical  approach 
is  scale  separation.  The  way  this  problem  has  been  addressed  in  the  literature  is  by  way  of 
dynamic  modeling,  which  provides  for  adaptive  selection  of  the  so-called  Smagorinsky  constant 
(see,  e.g.,  Germano  et  al.,  1991;  Moin  et  al.,  1991;  Piomelli,  1998). 

In  the  variational  multiscale  approach,  scale  separation  is  invoked  ab  initio.  A  space-time 
formulation  of  the  incompressible  Navier-Stokes  equations  is  employed.  From  the  discrete 
point  of  view,  this  leads  to  the  time-discontinuous  Galerkin  method  on  space-time  slabs.  This 
procedure  proves  convenient  and  obviates  the  need  to  consider  specific  time  discretization  (i.e., 
ODE)  algorithms.  However,  it  is  a  straightforward  matter  to  develop  traditional  semi-discrete 
approaches  with  similar  properties  to  the  ones  described  here.  The  appendix  to  this  section 
illustrates  this  fact. 

Two  simple  generalizations  of  the  Smagorinsky  eddy  viscosity  model  which  act  only 
on  small  scales  are  considered.  One  is  completely  desensitized  to  large-scale  behavior,  the 
other,  partially  desensitized.  Parameters  are  again  estimated  by  way  of  the  approach  due  to 
Lilly  (1966,1967,1988).  Because  scale  separation  has  been  invoked  from  the  outset,  a  constant- 
coefficient  model  in  the  present  approach  is  speculated  to  have  validity  for  a  greater  variety  of 
flows  than  the  classical  constant-coefficient  Smagorinsky  model.  The  approach  is  summarized 
by  contrasting  it  with  the  work  of  the  Temam  group  (see,  e.g.,  Dubois,  Jauberteau  and  Temam, 
1993,1998).  How  the  present  approach  addresses  criticisms  of  the  classical  LES/ Smagorinsky 
method  and  the  identification  of  some  other  useful  properties  are  described. 

6.1.  Incompressible  Navier-Stokes  equations 

Let  H  be  an  open,  connected,  bounded  subset  of  Rd ,  d  =  2  or  3,  with  piecewise  smooth 
boundary  T  =  dfl;  fi  represents  the  fixed  spatial  domain  of  our  problem.  The  time  interval  of 
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interest  is  denoted  ]0,T[,  T  >  0,  and  thus  the  space-time  domain  is  Q  =  x  ]0,T[;  its  lateral 
boundary  is  denoted  P  =  T  x  ]0,T[.  The  setup  is  illustrated  in  Figure  31. 


Figure  31.  Space-time  domain  for  the  initial/boundary- value  problem. 

The  initial/boundary-value  problem  consists  of  solving  the  following  equations  for  u  :  Q 
Rd ,  the  velocity,  and  p  :  Q  — >  ffi,  the  pressure  (divided  by  density), 


dxi 

—  +  V  •  (tt  <g>  u)  +  Vp  = 
ot 

vAxi  +  f 

in  Q 

(380) 

V  •  u  = 

0 

in  Q 

(381) 

u  = 

0 

on  P 

(382) 

«(0+)  = 

«(CT) 

on 

(383) 

where  f  :  Q  — ■>  is  the  given  body  force  (per  unit  volume);  v  is  the  kinematic  viscosity, 

assumed  positive  and  constant;  m(0_)  :  S!  — >  is  the  given  initial  velocity;  and  (g>  denotes 

the  tensor  product  (e.g.,  in  component  notation  [u  tg)  v] tJ  =  UiVj).  Equations  (380)-(383) 
are,  respectively,  the  linear  momentum  balance,  the  incompressibility  constraint,  the  no-slip 
boundary  condition  and  the  initial  condition. 

The  setup  for  a  space-time  formulation  is  recalled  once  again.  When  a  function  is  written 
with  only  one  argument,  it  is  assumed  to  refer  to  time.  For  example,  u(t)  =  where  the 

spatial  argument  x  €  is  suppressed  for  simplicity.  Furthermore, 

u  (i^)  =  lim  u  (t  i  e)  Vi  G  [0,T]  .  (384) 

This  notation  allows  us  to  distinguish  between  m(0+)  and  m(0_),  the  solution  and  its  given 
initial  value,  respectively.  In  the  variational  formulation  of  the  initial/boundary- value  problem 
(see  Section  6.4),  (383)  will  only  be  satisfied  in  a  weak  sense.  The  notation  of  (383)  and  (384) 
is  also  conducive  to  the  generalization  of  the  formulation  to  the  discrete  case  in  which  the 
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numerical  problem  is  posed  in  terms  of  a  sequence  of  “space-time  slabs,”  where  the  solution 
may  be  discontinuous  across  the  slab  interfaces. 

For  mathematical  results  of  existence,  uniqueness  and  regularity,  see  Temam  (1984), 
Quarteroni  and  Valli  (1994),  and  references  therein. 

Solutions  of  (380)-(383),  for  the  case  in  which  v  is  very  small,  typically  give  rise 
to  turbulence,  an  inherently  chaotic  and  unpredictable  phenomenon.  Nevertheless,  some 
statistical  quantities,  such  as  particular  spatial  and  temporal  averages,  are  deterministic  and, 
in  principle,  computable. 

6.2.  Large  Eddy  Simulation  (LES) 

The  unpredictability  of  turbulence  suggests  reformulating  the  initial/boundary-value  problem 
in  terms  of  averaged  quantities.  In  Large  Eddy  Simulation  (LES)  a  spatial  averaging  procedure 
is  employed.  For  example,  let 

g(x,y)u(y,t)  dy  (385) 

-Da  (X) 

in  which 

g{x,  y )  =  g(x  -  y)  (homogeneity)  (386) 

and 

1=  J  g(x,y)  dy  (387) 

Da(X) 

where  u  is  the  filtered  velocity  and  g  is  the  filter  having  support  in  DA{x)  C  il.  a  neighborhood 
of  x  C  n.  (See  Fig.  32  for  a  schematic  of  a  candidate  filter.) 


x_ 

A 


Figure  32.  Typical  filter  function  for  LES. 

The  size  of  DA(x)  is  characterized  by  A,  the  filter  width.  There  are  various  possibilities  for 
DA(x).  For  example,  a  possible  difinition  is 

DA(x)={yCRd  \p(x,y)<  A/2}  .  (388) 
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Figure  33.  The  effect  of  filtering. 


The  distance  function,  p,  may  be  defined  in  terms  of  the  Euclidean  norm,  in  which  case  D/\{x) 
is  an  open  ball  of  radius  A/2  centered  at  x. 

The  effect  of  filtering  is  schematically  illustrated  in  Figure  33. 

The  filtered  field,  U,  is  commonly  referred  to  as  the  large,  coarse,  or  resolvable  scales.  It  is 
assumed  adequate  to  represent  the  larger  structures  of  the  flow.  The  difference  between  u  and 
U,  that  is, 


u'  =  u  —  u  (389) 

is  the  rapidly  fluctuating  part  of  u ;  u'  is  commonly  referred  to  as  the  small,  fine,  or  unresolvable 
scales.  Although  primary  interest  is  in  computing  It,  due  to  the  nonlinear  nature  of  the  Navier- 
Stokes  equations,  the  effect  of  u'  on  U  cannot  be  ignored. 

Remark 

The  homogeneous  structure  of  the  filter  results  in  the  commutativity  of  spatial 
differentiation  and  filtering,  a  property  exploited  in  the  derivation  of  the  equations  governing 
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filtered  quantities.  However,  in  order  to  obtain  the  filtered  equations  corresponding  to  (380)- 
(383),  the  filtering  operation  needs  to  be  performed  for  all  ce  G  fi,  and  in  particular  for 
x  G  0\ClA,  where 

=  {x  |  Da(x)  c  H}  .  (390) 

In  this  case,  the  support  of  the  filtering  operation  extends  beyond  the  boundary  of  H.  For 
a  schematic  illustration,  see  Figure  34.  Clearly,  this  creates  mathematical  ambiguities.  (Note 
that  this  problem  does  not  arise  in  cases  of  domains  with  periodic  boundary  conditions.) 
An  alternative  approach  is  to  reduce  the  size  of  the  filter  as  the  boundary  is  approached,  but 
this  too  creates  mathematical  complications.  This  issue  will  not  be  addressed  further  herein.  It 
needs  to  be  emphasized,  though,  that  it  is  an  important  issue.  (For  recent  literature  concerning 
this  problem,  see  Galdi  and  Layton,  2000;  John  and  Layton,  1998;  Layton,  1996;  Ghosal  and 
Moin,  1995.) 

6.2.1.  Filtered  Navier-Stokes  equations 


du 

—  +  V  •  (u  0  u)  +  Vp  = 

vAu  +  f 

in  Q 

(391) 

V  •  u  = 

0 

in  Q 

(392) 

u  = 

0 

on  P 

(393) 

m(0+)  = 

m(0") 

on  fi  . 

(394) 

The  nonlinear  term  in  (391)  gives  rise  to  a  closure  problem :  how  to  compute  u®ul  This 
necessarily  entails  some  form  of  approximation.  To  this  end,  define  the  subgrid-scale  stress 

T  =  u<g>u  —  u<S>u  .  (395) 

In  terms  of  the  subgrid-scale  stress,  the  filtered  momentum  equation  (391)  is  rewritten  as 

Q'ui  _ 

— — I-  V  '(u®u)  +  Vp  =  vA u  +  V  ■  T  +  /  .  (396) 

at 

Thus,  T  needs  to  be  modeled  to  close  the  system.  To  be  more  precise,  only  the  deviatoric  part 
of  T,  namely, 

dev  T=  T-  (|  tr  T)  7  ,  (397) 

where  I  is  the  identity  tensor,  needs  to  be  modeled,  and  the  dilatational  part,  |  tr  T,  may  be 
subsumed  by  p. 

6. 3.  Smagorinsky  closure 

The  classical  and  most  widely  used  closures  are  based  on  the  Smagorinsky  eddy  viscosity 
model  Smagorinsky  (1963): 

Ts  =  2  uTVsu  (398) 

where 

vT  =  (CsA)2  |Vsu|  (399) 

Vsu  =  1(V«+(Vm)t)  (400) 

|VS«|  =  (2  Vsu-Vsu)1/2  (401) 
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and  Cs  is  referred  to  as  the  Smagorinsky  constant.  Note,  Ts  is  deviatoric,  i.e.,  Ts  = 
dev  Ts- 

Various  criticisms  have  been  lodged  against  the  Smagorinsky  model  (see,  e.g.,  Germano  et 
al.,  1991;  Piomelli,  1998).  Typical  of  these  are: 

1.  Ts  does  not  replicate  the  asymptotic  behavior  of  T  near  walls,  in  particular,  T s  does 
not  vanish  at  walls. 

2.  Values  of  Cs  obtained  from  the  decay  of  homogeneous  isotropic  turbulence  tend  to  be 
too  large  for  other  situations,  such  as  in  the  presence  of  mean  shear. 

3.  Ts  produces  excessive  damping  of  resolved  structures  in  transition,  resulting  in 
incorrect  growth  rate  of  perturbations. 

As  a  result  of  these  shortcomings,  many  modifications  have  been  proposed,  such  as  wall 
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functions,  intermittency  functions,  etc.  Perhaps  the  most  notable  achievement  is  the  dynamic 
subgrid-scale  model  (Germano  et  al.,  1991),  in  which  it  is  assumed  that  Cs  is  a  function  of 
space  and  time,  that  is 


Cs=Cs(x,t)  .  (402) 

The  identification  of  Cs  is  performed  adaptively  by  sampling  the  smallest  resolved  scales  and 
using  this  information  to  model  the  subgrid  scales.  The  dynamic  model  has  been  applied  to  a 
variety  of  flows  and  improved  results  have  been  obtained  in  most  cases.  For  a  recent  review  of 
the  state-of-the-art  and  assessment,  see  Piomelli  (1998).  It  is  a  widely  held  opinion  that  any 
proposal  of  a  new  LES  model  based  on  the  Smagorinsky  concept  must  address,  at  the  very 
least,  the  shortcomings  delineated  above. 

6.3.1.  Estimation  of  parameters  The  Smagorinsky  parameters  Cs  and  A  may  be  determined 
by  a  procedure  due  to  Lilly  (1966,1967,1988).  In  Lilly’s  analysis  it  is  assumed  that  turbulent 
kinetic  energy  dissipation  and  dissipation  produced  by  the  Smagorinsky  model  are  in  balance. 
The  limit  of  resolution  is  assumed  to  fall  in  the  Kolmogorov  inertial  subrange  and  |Vsu| 
is  determined  by  spectral  integration.  This  enables  quantification  of  Cs  A  and  vt-  A  brief 
summary  of  the  steps  involved  follows. 


In  E{k) 


Figure  35.  Kolmogorov  energy  spectrum. 

Consider  Figure  35.  E(k)  is  the  spectral  amplitude  of  kinetic  energy,  defined  as  the  integral 
over  surfaces  of  spheres  in  wave-number  space  parametrized  by  the  radius  k.  In  the  inertial 
subrange 

E(k)  =  a£2/3k~5/3  (403) 
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where  a  is  the  Kolmogorov  constant  and  £  is  the  turbulent  dissipation.  |Vstt|  may  be 
determined  from  the  relation 


k 

k2E(k)  dk 


(404) 


where  k  corresponds  to  the  resolution  limit,  which  is  the  cut-off  wave  number  for  spectral 
discretization.  Equating  turbulent  kinetic  energy  dissipation  with  dissipation  produced  by  the 
Smagorinsky  model,  and  evaluating  (404)  using  (403),  yields 


£ 


Ts  ■  Vsu 


2(CSA)2  |Vsm|  (Vsu-Vsu) 
(CSA)2  |VS«|3 


3/2 

k2e 


from  which  it  follows  that 


and 


(405) 


(406) 


(407) 


Assuming  k  =  const  h  1,  where  h  is  the  mesh  parameter,  it  follows  from  (406)  and  (407) 
that 


Cs  A  =  O(h)  (408) 

and 

vT  =  0(h4/3)  .  (409) 


Remarks 

1.  Assuming  A  =  h  =  nk"1  and  a  =  1.4,  it  follows  from  (406)  that  Cs  —  0.18.  However, 
smaller  values,  often  around  0.10,  are  often  used  in  practice  (Germano  et  al.,  1991  ). 

2.  The  excessive  damping  of  resolved  structures  may  be  explained  by  (398)  and  (409).  An 
0(h4/3)  viscosity  acts  on  all  scales  present.  It  is  known  from  the  analysis  of  artificial 
viscosity  methods  that,  even  for  linear  model  problems,  an  0(h 4/3)  artificial  viscosity 
results  in  convergence  that  is  at  most  0(h 4/3)  in  L2  and  0(h2^3)  in  H1.  This  is  indeed 
very  slow  convergence  and  is  deemed  by  most  analysts  as  unacceptable.  Furthermore, 
these  results  are  probably  optimistic  for  the  (nonlinear)  Navier-Stokes  equations.  The 
physical  design  of  the  Smagorinsky  model  results  in  the  correct  extraction  of  total 
kinetic  energy,  but  the  flaw  seems  to  be  that  the  extraction  of  kinetic  energy  occurs  in 
all  scales  and,  in  particular,  in  the  so-called  “resolved  scales.” 
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3.  The  analysis  described  assumed  an  isotropic  discretization.  The  anisotropic  case  has  also 
been  addressed  by  Lilly  (1988).  See  also  Scotti  and  Meneveau  (1993),  Scotti,  Meneveau 
and  Fatica  (1997). 

4.  The  concepts  of  filtering  and  the  filtered  equations  involve  subtleties  not  described  here. 
For  more  extensive  discussion  of  these  issues,  the  interested  reader  may  consult  Carati, 
Winckelmans  and  Jeanmart  (2001),  Winckelmans  et  al.  (2001),  and  Winckelmans, 
Jeanmart  and  Carati  (2002). 

5.  It  has  been  observed  that  some  common  filters  are  isomorphisms  (e.g.,  the  Gaussian 
filter).  In  these  cases,  u  can  be  reconstructed  from  u.  Consequently,  the  filtered 
equations  contain  the  same  information  as  the  Navier-Stokes  equations.  A  formulation 
of  LES  which  circumvents  the  use  of  the  filtered  equations,  and  associated  conceptual 
difficulties,  is  presented  next. 

6.f.  Variational  multiscale  method 

6.4.I.  Space-time  formulation  of  the  incompressible  Navier-Stokes  equations  Consider  a 
Galerkin  space-time  formulation  with  weakly  imposed  initial  condition.  Let  V  =  V(Q)  denote 
the  trial  solution  and  weighting  function  spaces,  which  are  assumed  to  be  identical.  Assume 
U  =  {u,p}  G  V  implies  u  =  0  on  P  and  fnp(t)  dfl  =  0  for  all  t  €  ]0,T[.  The  variational 
formulation  is  stated  as  follows: 

Find  U  €  V  such  that  VW  =  {to,  q}  €  V 

B(W,U)  =  (W,F)  (410) 


where 

B(W,U)  =  (w(T-),u(T~))n  -  q  ~  (VtiMi®u)Q 

+  (q,  V  •  u)q  -  (V  •  w,p)q  +  (Vsw,2 vVsu)Q  (411) 

and 

(W,  F)  =  (w,  f)Q  +  (w(0+),u(0~))n.  (412) 

This  formulation  implies  weak  satisfaction  of  the  momentum  equations  and  incompressibility 
constraint,  in  addition  to  the  initial  condition.  The  boundary  condition  is  built  into  the 
definition  of  V. 

Remarks 

1.  u(0_)  is  viewed  as  known  when  computing  the  solution  in  Q. 

2.  The  standard  weak  formulation  corresponding  to  the  discontinuous  Galerkin  method 
with  respect  to  time  is  obtained  by  replacing  [0,T]  by  [tnitn+ 1],  n  =  0,1,2,...  and 
summing  over  the  space-time  slabs 

Qn  =  ilx]tn,tn+1[  .  (413) 

In  this  case,  (410)-(412)  are  viewed  as  the  variational  equations  for  a  typical  slab. 
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3.  The  conditions  V  •  u  =  0  on  Q  and  u  =  0  on  P  imply 

(Vw,  u  ®  m)q  =  0  .  (414) 

In  the  discrete  case  this  term  may  need  to  be  altered  to  preserve  this  property.  See 
Quarteroni  and  Valli  (1994),  p.435. 

Kinetic  energy  decay  inequality : 

Substitution  of  U  for  W  in  (410)  leads  to  the  inequality 

IIKnlln  +  2^  ||  Vs'u||g  <  §||t*(0")||*  +  (u,f)Q 
from  which  follows 

5ll«(T-)||n  +  Hivs«||^  <  |||«(<r)||’  +  ^WfWl 

where  Cq  is  the  constant  in  the  Poincare  inequality: 

IMIn  <  Cn  ||Vsm||^  . 

6-4-2.  Separation  of  scales  Let 

V  =  V©V'  (418) 

The  reader  is  reminded  that  V  is  identified  with  a  standard  finite  element  space  and  V' 
is  oo-dimensional.  In  the  discrete  case,  V'  can  be  replaced  with  various  finite-dimensional 
approximations,  such  as  hierarchical  p-refinement,  bubbles,  etc.  In  any  case,  (418)  enables 
decomposition  of  (410)  into  two  sub-problems: 

B(W,U+U')  =  ( W,F ) 

B(W',U  +  U')  =  C W',F ) 

where 

u  =  u  +  u' 

w  =  W+W' 

in  which  U,  W  G  V  and  U' ,  W'  G  V'. 

It  is  important  to  realize  that  although  the  use  of  the  over-bar  and  prime  notations  continues 
to  connote  large  and  small  scales,  the  meaning  is  quite  different  than  for  the  classical  LES 
formulation  considered  previously.  Here  U  and  U'  may  be  thought  of  as  “projections”  of  U 
onto  V  and  V',  respectively.  The  terminology  “projections”  is  used  loosely  because  U  and  U' 
are  obtained  from  U  by  solving  coupled  nonlinear  problems,  viz., 

B(W,U  +  U')  =  B(W,U)  (423) 

B{W\U  +  U')  =  B(W',U)  .  (424) 

Consequently,  it  is  not  possible  to  identify  a  simple  filtering  relationship  between  U  and  U, 
such  as  (385)- (386).  Nevertheless,  U  represents  the  part  of  U  which  lives  in  V,  and  thus 
clearly  is  a  large-scale  representation  of  U.  Likewise,  U'  is  a  small-scale  representation  of  U. 
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The  relationship  between  the  present  U  and  its  filtered  counterpart  is  a  complex  mathematical 
problem.  Reference  may  be  made  to  Galdi  and  Layton  (2000)  and  John  and  Layton  (1998)  for 
initiatory  attempts  at  its  resolution. 

Let 

BAW^U.U')  =  4-B(W,U  +  eU')  (425) 

de  e=o 

—  d2  — 

B2(W  ,U  ,U')  =  —B{W  ,U  +  eU')  .  (426) 

£— 0 

With  these,  write 

B(W,U  +  U')  =  B{W,U)  +  Bi(W,U,  U')  +  ±B2{W,U,U')  (427) 

where 

±B2{W,U,U')  =  -(Vi»,u'®u')q  (428) 

and 

B1(W,U1  U')  =  (w{T-),u'{T-))^-{^,u^j  —  (Vw,u®  u'  +  u'  ®u)Q 

+  (9,V-w,)Q-(V-tn,p,)Q  +  (Vs^, 2i/W)0  .  (429) 

Bi(W,  U,  U')  is  the  linearized  Navier-Stokes  operator.  With  the  aid  of  (427)-(429),  rewrite 
(419)  and  (420)  as 

S(W,C7)  +  Ri(W,C7,i7')  =  (VuT,  u'  (g)  u')Q  +  (W\  F)  (430) 

B^WW^')-  {Vw',u'  ®u')Q  =  —  [B(W',U)  —  (W' ,F)]  .  (431) 

This  amounts  to  a  pair  of  coupled,  nonlinear  variational  equations.  Given  the  small  scales  (i.e., 
U1),  (430)  enables  solution  for  the  large  scales  (i.e.,  U).  Likewise,  the  large  scales  drive  the 
small  scales  through  (431).  Note,  the  right  hand  side  of  (431)  is  the  residual  of  the  large 
scales  projected  onto  V'. 

Remark 

The  subgrid-scale  stress,  T  (see  (395)),  may  be  decomposed  into  the  Reynolds  stress,  cross 
stress  and  Leonard  stress.  In  the  variational  formulation,  the  analogs  of  these  terms  are: 

(VuT,  u'  <S>  u') q  (Reynolds  stress)  (432) 

(VuJ,U(8)  v!  +  v!  ®  u)q  (Cross  stress)  (433) 

CVw,u(E)u)q  —  (Vw,u®u)q  =  —  (Vw',u®u)q  (Leonard  stress)  (434) 

Note  that  (432)- (433)  appear  in  (430),  and  (434)  appears  in  the  right-hand  side  of  (431).  Up 
to  this  point  our  results  are  exact ,  that  is,  nothing  has  been  omitted  and  no  modeling  has 
been  performed. 
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6.4-3.  Modeling  of  subgrid  scales  Observe  that  closure  is  not  the  motivation  for  modeling 
in  the  case  of  the  variational  multiscale  formulation.  In  fact,  there  is  no  closure  problem 
whatsoever.  The  need  for  modeling  is  due  simply  to  the  inability  of  typical  discrete 
approximations  to  properly  represent  all  necessary  scales.  This  is  viewed  as  a  conceptual 
advantage  of  projection/variational  methods  over  the  classical  filtered  equation  approach. 

Thus,  add  (Vw',R's)q  to  (431),  where 

R's  =  2 v'TVsu’  .  (435) 

(Specification  of  u'T  will  be  postponed  for  the  moment.)  The  end  result  is 

B'(W',U,  U')  =  -  [ B{W',U )  -  (W\  F)]  (436) 

where 

F'(Ty',t7,t/')  =  Fi(Ty',t7,[/')  -  (Vtu'^'  ®u')Q  +  (ysw',2v'TVsu')Q  .  (437) 

This  is  the  modeled  small-scale  equation  which  replaces  (431).  The  large-scale  equation,  (430), 
remains  the  same,  that  is,  there  is  no  modeling. 

Remark 

Collis  (2001)  has  presented  an  interpretation  of  this  formulation  that  is  more  consistent 
with  traditional  turbulence  modeling  concepts.  He  assumes 

u  =  u  +  u'  +  u" 
w  =  w  +  w'  +  w" 

and 

V  =  V  ©  V'  ©  V"  (440) 

where  V  and  V'  are  the  same  finite-dimensional  spaces  considered  here  and  V"  is  infinite¬ 
dimensional.  Thus  the  discrete  approximation  amounts  to  simply  omitting  V".  The  idea  is 
schematically  illustrated  in  Figure  36.  Collis  (2001)  may  also  be  consulted  for  clarifications  to 
the  formulation. 

6.4.4.  Eddy  viscosity  models  Two  candidate  definitions  of  v'T  will  be  considered.  In  both 
cases,  parameters  will  be  estimated  by  way  of  the  Lilly  analysis  (see  Section  6.3.1).  However, 
this  time  there  are  two  relevant  wave-number  scales:  k,  the  resolution  limit  of  the  space  V, 
and  k' ,  the  resolution  limit  of  the  space  V  =  V  ©  V'.  The  interpretation  of  k  is  assumed  to  be 
similar  to  before  (see  Section  6.3.1).  Figure  37  schematically  contrasts  the  present  situation 
with  classical  LES. 

Note  that  in  generalizing  the  Lilly  analysis  to  the  current  situation,  it  is  necessary  to  calculate 
spectral  integrals  over  the  interval  [k,  k']  (see  Fig.  38).  Consequently,  it  is  assumed  this  interval 
lies  entirely  within  the  inertial  subrange. 

The  assumed  forms  of  v'T  are: 

4  =  (CgA1)2  |Vsw'|  (441) 
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Figure  36.  Collis  interpretation  of  the  variational  multiscale  formulation. 
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Figure  37.  Above:  in  classical  LES,  the  subgrid-scale  stress,  Ts,  acts  on  all  scales  present,  namely 
k  £  [0,  k].  Below:  in  the  multiscale  model,  the  model  acts  only  on  the  small  scales,  namely  k  £  [fc,  k']. 


and 


4  =  (C'^A')2|Vsm|  .  (442) 

In  the  first  case,  4  depends  exclusively  on  small-scale  velocity  components.  In  the  second 
case,  4  depends  on  the  large-scale  components. 

As  before,  |Vsw|  is  evaluated  through  (403)  and  (404);  the  small-scale  counterpart  is 
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Figure  38.  Kolmogorov  energy  spectrum.  The  interval  of  small  scales  is  assumed  to  lie  within  the 

inertial  subrange. 


evaluated  with  (403)  and 


k' 

k2E(k)  dk  . 


With  this  result  in  hand,  calculations  similar  to  (405)  yield: 


/  9  \  3/4 

C'gA'  =  [£-)  k-1  [(fc7fc)4/3-l]-3/4 

/  9  \  ^/^ 

C'gA'  =  k-1  [(fc7fc)4/3-l]-1/2 

\3ay 


(443) 


(444) 

(445) 


corresponding  to  (441)  and  (442),  respectively.  For  a  given  discretization,  from  which  k  and  k1 
can  be  determined,  (444)  and  (445)  provide  estimates  of  C’gA'  for  the  two  cases  considered. 


Remarks 

1.  Note  that  fixing  k  in  (444)  and  (445)  implies  that 

k'  T  I  ^  C'gA1  I  T  (446) 

for  both  cases,  confirming  the  intuitively  obvious  result  that  the  inclusion  of  more  small 
scales  reduces  the  size  of  C's  A',  and  vice  versa.  If  A'  =  7 rfc-1,  as  in  Remark  1  of  Section 
6.3.1,  (444)  and  (445)  show  that  C’s  is  just  a  function  of  the  Kolmogorov  constant,  a, 
and  the  ratio  fc'/fc. 
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2.  The  argument  can  also  be  reversed:  For  fixed  k  and  C's A'  =  Cs A  (i.e. ,  the  value  for 
the  classical  Smagorinsky  model),  what  is  k'l  The  answer  to  this  question  estimates  the 
size  of  V'  in  comparison  with  V.  It  turns  out  that  for  both  (444)  and  (445), 

k'  =  23/ 4 k  «  1.7 k  .  (447) 

Roughly  speaking,  the  resolution  limit  needs  to  be  almost  twice  as  great.  If  it  is  assumed 
that 

k'/k  =  h./ti,  (448) 

(447)  implies 

h'  =  ft/ 23/4  «  0.6ft  ,  (449) 

which  amounts  to  a  finer  mesh  by  almost  a  factor  of  two  in  each  direction. 

3.  An  even  more  attractive  option  may  be  the  analytical  determination  of  U'  and  its  a 
priori  elimination  from  the  large-scale  equation,  that  is,  (430).  This  gives  rise  to  a 
nonlinear  stabilized  method.  Note  that  for  (447)- (449),  as  before, 

CSA'  =  0(h) 

4  =  0(ft4/3) 

However,  keep  in  mind  that  R's  only  acts  on  the  small  scales. 

6-4-5.  Precis  of  results  For  convenience,  the  main  LES/ variational  multiscale  equations  are 
summarized  here: 

Large  scales 

B(W,  U)  +  Bi  (W,  U,  U')  =  (Vw,u'  ®u')Q  +  (W,F)  (452) 

Small  scales 

Bi(W',U,  U')  —  (Vu/,  u'  <g>  u')Q  =  ~[B(W',U)-  (W',F)]  (453) 

Modeled  small  scales 

B\W',U,  U')  =  -  [B(W\  U)  -  (W\  F)]  (454) 

where 

B'(W\U,U')  =  Bi(W',U,  U')  -  (Vw\u'  ®u’)Q  + (Vs w',2v'tVsu')q  (455) 

and 

4  =  (C's A')2  |  Vs«'|  (456) 

or 

4  =  (C^A')2|Vsm|  .  (457) 


(450) 

(451) 
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Modeled  system 

A  concise  way  of  writing  the  combined  system  of  (452)  and  (454)  is 

B(W,U)  =  ( W,F ) 

where 

B(W,  U)  =  B(W,  U)  +  (Vs id',  2z4  W)q 

Kinetic  energy  decay  inequality  for  the  modeled  system 

This  follows  immediately  from  (458)  with  W  replaced  by  U : 

i||M(T-)||J  +  2H|Vs«||!+  (2 4)1/2W  2  <  I||m(0-)||^  +  (m,/)q  (460) 

which  leads  to 

i||tx(T-)||J!  +  ^||VsW||^+||(24)1/2w||2Q<  |||«(0-)||i  +  ^II/I|q  .  (461) 

Remarks 

1.  Note,  if  U  is  an  exact  solution  of  the  Navier-Stokes  equations,  then 

B(W,U)  =  (W,F)  VW  e  V  .  (462) 

In  particular, 

B(W',U)  =  (W',F)  VW'  £  V'  .  (463) 

Consequently,  for  both  the  exact  and  the  modeled  small-scale  equations  (i.e.,  (453)  and 
(454),  respectively),  it  follows  that  U'  =  0.  From  (462),  it  also  follows  that 

B(W,U)  =  (W,  F)  VW  e  V  •  (464) 

These  results  verify  that  (452)  is  identically  satisfied.  This  means  the  equation  governing 
the  large  scales  retains  its  consistency ,  or  residual  structure,  despite  the  modeling  of 
small  scales,  in  contrast  to  the  classical  LES  case  described  in  Section  6.2. 

2.  By  virtue  of  the  fact  that  U'  =  {u',p'}  €  V'  C  V,  u!  =  0  on  P.  Furthermore,  the 
small-scales  equations  (either  (453)  or  (454))  imply  V  •  u'  =  0  on  Q.  From  this  it  may 
be  concluded  that  u'  attains  the  correct  asymptotic  structure  near  walls  by  the  usual 
argument.  In  particular,  the  [1,2]  component  of  u'  (g>  u'  in  (452)  behaves  like  x\,  where 
X\  is  the  streamwise  direction  and  Xi  is  the  direction  normal  to  the  wall. 


(458) 

(459) 
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6.5.  Relationship  with  other  methods 

6.5.1.  Nonlinear  Galerkin  method  The  approach  has  both  similarities  to,  and  differences 
with,  the  work  of  Temam  and  colleagues.  The  similarities  are  that  both  approaches  invoke  scale 
separation  from  the  outset  and  both  employ  projected  forms  of  the  Navier-Stokes  equations, 
rather  than  the  filtered  versions  used  in  classical  LES.  The  salient  differences  are  as  follows: 
In  the  earlier  work  of  the  Temam  group  (see  Dubois,  Jauberteau  and  Temam,  1993  for  a 
representative  exposition) : 

(i)  In  the  equation  governing  large-scales,  the  classical  Reynolds  stress  term  is  neglected. 

(ii)  In  the  equation  governing  small  scales,  the  only  term  acting  on  small  scales  which  is 
retained  is  the  Stokes  operator  term.  In  particular,  the  cross-stress,  time-derivative  and 
small-scale  Reynolds  stress  terms  are  neglected. 

These  assumptions  are  strong  ones  and  can  be  shown  to  limit  the  applicability  of  the  approach 
to  discretizations  which  are  only  somewhat  coarser  than  required  for  a  fully-resolved  direct 
numerical  simulation  (DNS).  Otherwise  excessive  dissipation  is  encountered.  This  may  be 
contrasted  with  the  present  approach  in  which  all  these  terms  are  retained.  Obviously,  it  is 
assumed  here  to  be  important  to  retain  the  terms  omitted  in  the  Temam’s  group  earlier  work. 
In  addition,  a  model  is  added  to  the  equation  governing  the  small  scales  in  the  present  case. 

In  more  recent  work  (see  Dubois,  Jauberteau  and  Temam,  1998),  Temam  and  co¬ 
workers  retain  all  the  terms  omitted  in  their  previous  work  and  employ  an  adaptive  strategy 
in  space  and  time  to  resolve  small-scale  behavior.  This  seems  to  be  a  DNS  approach,  although 
Dubois,  Jauberteau  and  Temam  (1998)  characterize  it  as  somewhere  between  LES  and  DNS. 
The  present  approach  is  similar,  except  for  the  model  in  the  small-scale  equation.  It  is  felt 
appropriate  to  characterize  the  present  formulation  as  LES,  or  at  least  closer  to  LES  than  that 
of  Dubois,  Jauberteau  and  Temam  (1998). 

6.5.2.  Adaptive  wavelet  decomposition  Farge,  Schneider  and  Kevlahan  (1999)  present  a 
very  interesting  analysis  of  two-dimensional  turbulence  based  on  an  adaptive  wavelet 
decomposition  of  the  vorticity  into  coherent,  non-Gaussian  structures,  and  an  incoherent, 
Gaussian  background.  The  efficiency  of  the  wavelet  representation  is  illustrated  in  a  DNS 
computation  of  a  mixing  layer  in  which  the  time  evolution  of  the  coherent  part  only  involves 
8%  of  the  wavelet  coefficients.  The  remaining  coefficients,  associated  with  the  incoherent  part, 
are  simply  discarded  at  each  time  step.  In  order  to  maximize  data  compression,  the  velocity- 
vorticity  form  of  the  Navier-Stokes  equations  is  employed  in  preference  to  the  velocity-pressure 
form. 

The  wavelet  decomposition  at  each  time  step  is  an  example  of  a  priori  scale  separation,  as 
advocated  in  the  present  work.  However,  when  Farge,  Schneider  and  Kevlahan  (1999)  discuss 
modeling  they  do  so  in  the  context  of  the  filtered  form  of  the  equations,  as  in  classical  LES, 
and  are  faced  with  the  problem  of  representing  the  subgrid-scale  stress  T.  They  propose  the 
following  possibilities:  (i)  The  Smagorinsky  model;  (ii)  a  dynamic  generalization  in  which  the 
eddy  viscosity  vt  is  estimated  in  terms  of  the  enstrophy  fluxes  in  wavelet  space,  such  that  when 
energy  flows  from  large  to  small  scales  vt  will  be  positive  and  vice  versa  (i.e.,  backscatter); 
and  (Hi)  the  subgrid-scale  stress  T  modeled  as  a  Gaussian  forcing  term  proportional  to  the 
variance  of  the  incoherent  parts  of  the  vorticity  and  velocity.  In  modeling  T  in  these  ways  it 
seems  inevitable  that  at  least  some  of  the  shortcomings  associated  with  the  use  of  the  filtered 
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Perfectly-matched  layer 


Artificial  dissipation 


Figure  39.  The  perfectly-matched  layer. 


equations  will  eventually  be  encountered.  For  example,  T  may  act  too  strongly  on  the  coherent 
(i.e.,  resolved)  structures,  and  there  seems  no  a  priori  control  of  kinetic  energy  if  backscatter 
occurs.  However,  it  would  seem  possible  to  reformulate  the  numerical  procedure  in  terms  of 
the  variational  multiscale  method,  associating  coherent  structures  with  the  space  V  and  the 
incoherent  background  with  V'  ■  In  this  way  some  difficulties  might  be  circumvented,  but  the 
efficiency  of  the  method  might  also  be  compromised. 

6.5.3.  Perfectly-matched  layer  in  electromagnetics  and  acoustics  Consider  an  exterior  infinite 
domain  fi  C  Rd,  d  =  1,2,  or  3.  In  electromagnetic  wave  propagation  there  is  interest  in 
solving  exterior  infinite-domain  problems  with  finite  elements.  A  finite  region  is  introduced, 
discretized  by  finite  elements  and  surrounded  by  a  truncated  boundary.  The  situation  is  seen 
to  be  the  same  as  for  the  acoustics  problem  considered  previously  in  Section  2,  but  here  there 
is  no  attempt  to  solve  the  far-field  problem  exactly  and  severe  restrictions  are  not  placed 
on  the  shape  of  the  external  region.  For  example,  the  external  region  need  not  be  amenable 
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to  a  separation-of-variables  solution,  as  is  required  for  the  Dirchlet-to-Neumann  formulation 
considered  in  Section  2.  In  order  to  absorb  out-going  waves  and  to  prevent  spurious  reflections, 
a  finite  layer  is  introduced  in  which  the  properties  of  the  medium  are  artificially  altered  to 
be  dissipative.  The  dissipation  in  the  layer  is  designed  to  match  that  of  the  inner  region  at 
the  interface,  which  is  typically  non-dissipative,  and  to  progressively  increase  further  into  the 
layer.  To  ensure  a  smooth  transition  between  the  inner  region  and  the  layer,  the  dissipation 
coefficient  is  assumed  to  be  zero,  and  to  have  zero  derivative,  at  the  interface.  A  quadratic 
variation  is  seen  to  be  a  simple  way  to  achieve  the  desired  dissipative  properties  in  the  layer, 
and  is  one  often  used  in  practice  (see  Fig.  39.)  The  dissipative  layer  is  said  to  be  “perfectly 
matched.”  The  simplicity,  generality  and  effectiveness  of  the  approach  have  led  to  its  popularity 
in  electromagnetics,  despite  the  lack  of  a  rigorous  theory  to  support  the  precise  variation  and 
amplitude  of  artificial  dissipation.  Harari,  Turkel  and  Slavutin  (2000)  may  be  referred  to  for  a 
penetrating  analysis  and  references  to  important  literature.  The  perfectly-matched  layer  has 
features  in  common  with  ideas  and  techniques  used  in  turbulence  as  evidenced  by  the  following. 


Figure  40.  DNS  Results  for  homogeneous  turbulence,  Domaradzki  et  al.  ,  Re\  =  70  (Taylor  microscale 

Reynolds  number). 


Remarks 

1.  The  concept  of  “spectral  eddy  viscosity,”  introduced  by  Heisenberg  (1948),  has 
undergone  extensive  study  in  the  turbulence  literature  (see,  e.g.,  Kraichnan,  1976; 
Domaradzki,  Liu  and  Brachet,  1993;  Domaradzki  et  al.,  1994;  McComb  and  Young, 
1998).  The  idea  is  to  take  the  spatial  Fourier  transform  of  the  Navier-Stokes  equations, 
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Figure  41.  DNS  results  for  wall-bounded  flow,  Domaradzki  et  al.  ,  ReT  =  210  (Reynolds  number  based 

on  friction  velocity). 


and  eliminate  the  pressure  by  projection  onto  the  space  of  divergence-free  modes, 

—7^-  =  —  vk2ii(k)  +  N(k\u(q),  u(p))  (465) 

Here,  fc,  p ,  and  q  are  wave  vectors  and  N(k\ii(q),  u{p))  is  the  Fourier  transform  of 
the  nonlinear  term.  Equation  (465)  is  then  multiplied  by  the  complex  conjugate  of 
the  Fourier  coefficient,  namely  u*(k),  and  integrated  over  spherical  surfaces  of  radius 
k  =  |fc|  in  wave- vector  space,  to  derive  an  evolution  equation  for  the  energy  spectrum, 
E(k), 

dE}.k^  =  ~2vk2E{k)  +  T(k\u(q),u(p))  (466) 

at 

where  T(k\u(q),u(p))  is  referred  to  as  the  energy  transfer  term.  It  represents  the 
energy  transferred  from  the  spherical  surface  of  radius  fc,  to  Fourier  modes  with  wave 
numbers  p  and  q.  This  equation  can  be  given  a  numerical  analysis  interpretation  by 
introducing  a  cutoff,  /cc,  and  thinking  of  it  as  representing  the  limit  of  resolution  of 
a  discretization,  in  this  case,  spectral  truncation.  All  Fourier  modes  with  wave-vector 
amplitude  larger  than  kc  are  viewed  as  missing.  The  interesting  case  is  energy  transfer 
for  wave  numbers  k  <  kc.  It  makes  sense  to  split  the  energy  transfer  term  into  two  parts 
as  follows: 

dEg^  =  -2rk2E(k) +  T<(k\u{q),u{p)) +  T>(k\u{q)1ii{p))  (467) 

where  T<(k\ii(q),u(p))  assumes  p  and  q  are  smaller  in  magnitude  than  kc ,  hence  T<  is 
exactly  included  in  the  numerical  model,  and  T>(fc|u.(q),  u{p))  assumes  at  least  one  of  p 
and  q  is  larger  in  magnitude  than  kc,  hence  T>  is  omitted  in  the  numerical  model.  One 
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Figure  42.  Direct-interaction  approximation  (DIA)  results  for  homogeneous  turbulence  from 

Kraichnan. 


view  of  the  modeling  problem  in  turbulence  is  to  approximately  represent  this  omitted 
term.  Since  it  typically  represents  a  flow  of  energy  from  the  surface  of  radius  k  to  modes 
beyond  the  cutoff,  it  has  been  customary  to  rescale  it  as  a  spectral  eddy  viscosity, 
that  is, 


vt  (k)  = 


-r>(fc|M(q),w(p)) 

2  k2E{k) 


(468) 


Theoretical  studies  of  Kraichnan  (1976),  assuming  infinite  Reynolds  number,  suggested 
the  form  of  vr{k)  involves  a  plateau  at  low  wave  number  and  cusp  at  higher  wave 
number,  peaking  at  the  cutoff  (see  Fig.  42).  Initial  studies  of  DNS  data  bases  (well- 
resolved  turbulent  solutions  at  low  Reynolds  number)  confirmed  the  presence  of  the 
cusp,  but  the  plateau  was  at  zero,  rather  than  the  value  predicted  by  the  Kraichnan 
theory  (see  Fig’s.  40  and  41  from  Domaradzki,  Liu  and  Brachet,  1993;  Domaradzki  et 
al.,  1994).  The  issue  was  clarified  by  McComb  and  Young  (1998)  who  studied  a  higher 
Reynolds  number  flow  and  systematically  varied  the  location  of  the  cutoff  kc.  The  cusp 
is  clearly  apparent  in  all  cases,  but  the  plateau  also  emerges  at  a  very  low  value  of  the 
cutoff  (see  Fig.  43).  Keeping  in  mind  that  a  spectral  eddy  viscosity  model  acts  as  an 
artificial  viscosity  and,  due  to  the  predominance  of  the  cusp,  it  is  very  reminiscent  of 
the  cusp-like  structure  of  the  artificial  dissipation  employed  in  the  perfectly-matched 
layer.  Both  mechanisms  attempt  to  remove  information  which  would  be  transferred  to 
scales  present  in  the  numerical  approximations.  This  amounts  to  an  analogy  between 
wave-vector  space  behavior,  in  the  case  of  turbulence,  with  physical  space  behavior,  in 
the  case  of  the  perfectly-matched  layer.  The  latter  case  is  hyperbolic,  and  the  former 
represents,  in  aggregate,  an  almost  hyperbolic  phenomenon  in  wave-vector  space.  It  is 
interesting  that  similar  dissipative  mechanisms  are  used  to  represent  seemingly  very 
different  phenomena,  namely,  radiation  damping  in  the  electromagnetic  case,  a  linear 
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Figure  43.  DNS  results  from  McComb  &  Young  on  a  2563  mesli  at  Rex  =  190,  16.5  <  kc  <  112.5. 
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hyperbolic  mechanism  in  physical  space,  and  energy  transfers  due  to  nonlinearities  in 
turbulence,  which  become  nonlocal  interactions  in  wave- vector  space.  The  predominance 
of  the  cusp  in  both  cases  suggests  a  similar  mechanism  may  be  present.  In  fact,  it  can 
be  shown  in  the  case  of  turbulence  that  the  cusp  is  produced  by  linearized  interactions 
(i.e. ,  “cross  stress”  terms)  and  the  transfers  are  confined  to  a  spherical  layer  in  wave- 
vector  space  of  precisely  twice  the  radius  of  the  cutoff,  kc  (see  Hughes,  Wells  and  Wray, 
2003). 

2.  The  version  of  the  variational  multiscale  method  described  previously  may  be  thought  of 
as  an  approximation  of  spectral  eddy  viscosity  in  that  the  artificial  viscosity  inspired  by 
the  Smagorinsky  model  is  assumed  to  act  only  in  a  high  wave-number  layer  in  spectral 
space  (see  Fig.  44).  Likewise,  it  may  be  thought  of  as  a  wave- vector  interpretation  of  the 
perfectly-matched  layer,  although  no  attempt  has  so  far  been  made  to  attain  a  smooth 
match. 

3.  “Hyperviscosities”  in  spectral  formulations  (see,  e.g.,  Borue  and  Orszag,  1994,1996), 
which  are  also  cusp- like,  may  be  viewed  as  approximations  of  spectral  eddy  viscosity. 
As  pointed  out  by  Cerutti,  Meneveau  and  Knio  (2000),  however,  their  behavior  in 
physical  space  does  not  correlate  very  well  with  turbulent  energy  transfers.  When  inverse 
transformed  to  physical  space,  hyperviscosities  lead  to  higher-order  spatial  differential 
operators,  which  cannot  be  easily  implemented  in  many  numerical  methods  such  as,  for 
example,  finite  elements. 

6.5.^-  Dissipative  structural  dynamics  time  integrators  It  was  observed  many  years  ago  that 
it  was  important  to  introduce  dissipative  mechanisms  in  structural  dynamics  time  integrators 
to  prevent  spurious  oscillations  and  dangerous  “overshoot  phenomena”  in  high  frequencies, 
which  turn  out  to  be  equivalent  to  high  wave-vector  components  (see  Hughes,  1987,  Chapter 
9,  and  references  therein).  In  structural  dynamics  it  was  argued  that  time  integrators  could 
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Figure  44.  In  three-dimensional  turbulence,  energy  transfers  in  spectral  space  cascade  from  small  wave 
number  to  large  wave  number,  reminescent  of  wave  propagation  phenomena  in  exterior  problems  in 
physical  space.  In  the  variational  multiscale  method,  artificial  viscosity  is  incorporated  in  the  annular 
layer  [k,  k ']  to  approximate  the  energy-transfer  mechanism.  This  is  reminescent  of  the  perfectly- 

matched  layer. 
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At  IT 

Figure  45.  Spectral  radii  for  Newmark  methods  with  varying  /3. 


be  developed  to  remove  high-frequency  components  of  the  solution,  but  that  care  needed  to 
be  exercised  in  order  not  to  degrade  accuracy  in  low-frequency  components.  The  HHT  a- 
method  (see  Hilber,  Hughes  and  Taylor,  1977;  Hilber  and  Hughes,  1978)  and  Park’s  method 
(Park,  1975)  are  two  well-known  integrators  designed  to  satisfy  these  objectives.  For  further 
discussion  and  evaluation,  see  Hughes  (1987).  There  again  is  an  analogy  with  the  variational 
multiscale  method  in  that  dissipation  is  added  in  the  high  wave-number  regime  to  remove 
spurious  pile-up  of  energy  near  the  cutoff  (i.e. ,  limit  of  resolution)  whereas  it  is  avoided  in  the 
low  wave- number  regime  so  as  not  to  degrade  accuracy  in  the  well-resolved  modes.  There  are 
similarities  with  the  other  ideas  described  previously  mutatis  mutandis.  As  a  measure  of  the 
high-frequency  dissipation,  one  often  uses  the  spectral  radius  of  the  amplification  operator. 
A  value  of  1  indicates  no  dissipation  and  a  value  of  0  indicates  maximal  dissipation.  Figure 
45  illustrates  the  cusp-like  behavior  of  spectral  radii  of  some  dissipative  integrators  used  in 
structural  dynamics.  Note  the  absence  of  dissipation  in  low  frequencies  and  the  presence  of 
dissipation  in  high  frequencies. 

6. 6.  Summary 

In  this  section  a  variational  multiscale  approach  to  Large  Eddy  Simulation,  has  been  described. 
The  approach  may  be  contrasted  with  classical  LES  in  which  modeling  of  the  subgrid-scale 
stress  in  the  filtered  equations  is  performed  first,  and  scale  separation  is  accomplished  a 
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posteriori  by,  for  example,  the  dynamic  modeling  procedure. 

Criticisms  lodged  against  the  classical  LES/Smagorinsky  model  are  addressed  in  the  context 
of  the  present  approach.  Specifically,  even  with  a  positive,  constant-coefficient,  eddy  viscosity 
term,  the  present  approach: 

1.  Yields  the  correct  asymptotic  behavior  at  walls  for  the  classical  Reynolds  stress,  cross¬ 
stress,  Leonard  stress  and,  consequently,  the  subgrid-scale  stress.  In  particular,  no 
modeling  of  these  terms  is  performed. 

2.  Reduces  dissipation  in  the  presence  of  mean  shear. 

3.  Ameliorates  damping  of  resolved  structures  since  the  eddy  viscosity  term  acts  only  on 
fluctuations. 

Furthermore,  the  present  approach  de-emphasizes  the  role  of  the  Smagorinsky  constant 
in  producing  dissipation,  but  accentuates  dissipation  effects  associated  with  fluctuations. 
Consequently,  a  single  value  of  the  constant  is  anticipated  to  behave  in  a  more  satisfactory 
manner  for  a  variety  of  complex  flows. 

In  addition,  the  variational  equation  governing  large  scales,  which  is  unmodeled,  is 
identically  satisfied  by  all  exact  solutions  of  the  Navier-Stokes  equations.  This  is  often  referred 
to  as  the  “consistency  condition”  and  is  a  key  condition  for  obtaining  error  estimates.  This 
may  be  contrasted  with  the  classical  LES/constant-coefficient  Smagorinsky  model  case  which 
does  not  possess  this  property. 

On  the  other  hand,  the  eddy  viscosity  models  considered  herein  are  very  simple,  perhaps 
too  simple.  The  energy  transfers  from  low  wave  numbers  to  unresolved  modes,  present  in 
coarse  and  inviscid  LES  and  manifested  by  the  plateaus  in  Figures  42  and  43,  are  not  well 
represented  by  small-scale  viscous  mechanisms.  Consequently,  the  search  for  better  methods 
continues.  Nevertheless,  experience  so  far  with  the  present  method  has  been  remarkably  good. 
Results  with  spectral  formulations  of  homogeneous  isotropic  flows  and  turbulent  channel  flows 
were  presented  in  Hughes  et  al.  (2001),  Hughes,  Oberai  and  Mazzei  (2001)  and  Oberai  and 
Hughes  (2002).  Improved  results  with  a  dynamic  version  have  been  reported  on  in  Holmen 
et  al.  (2004),  and  Hughes,  Wells  and  Wray  (2003).  The  variational  multiscale  formulation 
of  LES  was  also  extended  by  Farhat  and  Koobus  (2002)  and  Koobus  and  Farhat  (2004)  to 
the  compressible  Navier-Stokes  equations  for  unstructured  finite  element  and  finite  volume 
discretizations.  Results  for  vortex  shedding  dominated  flows  are  presented  in  Koobus  and 
Farhat  (2004)  and  shown  to  match  experimental  data.  Other  works  which  present  results  for 
the  variational  multiscale  formulation  of  LES  and  variations  are  Winckelmans  and  Jeanmart 
(2001),  Jeanmart  and  Winckelmans  (2001),  Collis  (2002),  and  Ramakrishnan  and  Collis 
(2002, 2004a, 2004b, 2004c). 

A  summary  of  the  main  tenets  of  the  approach  described  herein  is  presented  as  follows: 

(i)  Variational  projection  in  preference  to  filtering.  This  obviates  the  closure  problem  and 
complex  issues  associated  with  filtering  and  inhomogeneous  wall-bounded  flows  are  also 
eliminated. 

(ii)  A  priori  scale  separation  in  preference  to  a  posteriori  scale  separation.  This  facilitates 
modeling  restricted  to  unresolved,  high-wave  number  phenomena  rather  than  all  wave 
numbers  as  in  classical  LES. 
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( in)  Modeling  confined  to  the  small-scale  equation  in  preference  to  modeling  within  the  large- 
scale  equation.  This  means  that  the  large-scale  equation  is  unmodeled  and  is  consistent 
in  the  weighted  residual  sense,  in  contrast  to  the  modeled  filtered  equations. 

6.7.  Appendix:  Semi- discrete  formulation 

In  this  appendix,  the  ideas  are  illustrated  in  terms  of  a  more  traditional  semi-discrete 
formulation  in  which  a  midpoint  rule  algorithm  is  used  for  time  discretization.  Let  un  and 
pn  represent  the  algorithmic  solution  at  time  tn.  The  time  step  is  denoted  At  =  tn+\  —  tn.  It 
proves  convenient  to  employ  the  jump  and  mean  value  operators,  viz., 

[«]  =  un+ 1  -  u„  (469) 

(it)  =  \  {un+ 1  +  un)  .  (470) 

The  algorithm  is  an  approximation  to  the  semi-discrete  variational  formulation,  given  as 
follows: 

(  du  \ 

—  J  -  ( Vw,u®u)n  +  (q,  V  •  u)n  -  (V  •  w,p)a  +  (Vsw,2vVsu)n  =  ( w,f)n  (471) 

Midpoint  rule 

-  (Vin,  (u)  ®  (u))n  +  (q,  V  •  (tx))n  -  (V  •  w,  {p))n 

+  (Vsw,2vVs(u))n  =  (w,(f))Q  (472) 

The  kinetic  energy  evolution  law  immediately  follows  from  (472)  by  replacing  w  with  (it)  and 
q  with  (p),  i.e. , 

\  llM"+illn  +  2vAt  llvs(M)llo  =  ^IMIn  +  ((“),(/)) a  ■  (473) 

In  addition, 

!ll«n+1||n  +  ^At  ||Vs(m)||q  <  3  ll^nlln  +  ^||(/)||n  (474) 

where  Cq  is  defined  by  (417). 

The  multiscale  procedure  is  developed  in  analogous  fashion  to  the  space-time  case  (see 
Section  6.4).  Let 

un  =  un  +  u'n  (475) 

Pn  =  Pn+p'n  (476) 

w  =  w  +  w'  (477) 

q  =  q  +  q '  (478) 

These  are  substituted  into  (472)  resulting  in  equations  governing  large  and  small  scales,  as 
before.  The  details,  which  are  straightforward,  are  left  to  the  interested  reader.  Modeling  also 
follows  the  ideas  developed  previously,  namely,  to  the  left-hand  side  of  (472)  add  the  term 

(VW,2 4V>'»q  ■  (479) 
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This  amounts  to  adding  viscous  dissipation  to  the  small  scales  equation.  The  modification  to 
the  kinetic  energy  identity,  (473),  is 


\  ll^n+illn  +  2^At  ||Vs(m)||q  +  At  |  |(2^)1,/2Vs(m') 
from  which  follows: 

\  ||«„+1||n  +  vAt  ||Vs(m)||^  +  At  I  \{2u'Tf/2Vs(u') 


£  +  Ai«t*),</»n  (480) 


<5 


u. 


2 

■nllfl  ' 


^ll</>llS, 


(481) 


Remark 

As  remarked  previously,  the  term  (V  (w),  ( u )  ®  (tt))n  may  need  to  be  altered  in  the  discrete 
case  in  order  to  achieve 


(V(u),(u)®(«»n=0  . 
See  Quarteroni  and  Valli  (1994),  p.435. 


(482) 
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